Skip to content

optim and singularity

5 messages · Tom Purucker, Berend Hasselman, Morway, Eric

Morway, Eric
# 
Hello,

I was unable to find clues to my problem in ?optim.  Using the data and code
below, I get an error ("system is exactly singular") when a particular line
of code is left in, but have found that 'optim' works when I comment it out. 
The line of code in question is after the closeAllConnections() line of code
and contains a call to "na.approx" from the zoo package.  What's confusing
to me is why optim works without it versus with it?  Any clues would be very
helpful in moving forward.  The overall goal of the code below is to
minimize the residuals between the blue circles and green squares by
adjusting the limits of the secondary y-axis.

Thanks, Eric

library(zoo)
temp.dat<-read.table(textConnection("Well Meas_Date WTD ECgw GW_Project_Area
Region Avg.EM.Survey.Value Avg.Soil.Paste.EC SD_EM_Survey
1 4/12/1999 1.75 2.27 LOWER-UPSTREAM US NA NA NA
1 5/11/1999 1.24 5.04 LOWER-UPSTREAM US NA NA NA
1 5/27/1999 1.27 4.45 LOWER-UPSTREAM US NA NA NA
1 6/3/1999 1.27 4.09 LOWER-UPSTREAM US 3.347069 3.126667 0.6347013
1 6/11/1999 1.52 2.84 LOWER-UPSTREAM US NA NA NA
1 6/18/1999 1.19 2.34 LOWER-UPSTREAM US NA NA NA
1 6/24/1999 1.21 2.48 LOWER-UPSTREAM US NA NA NA
1 6/30/1999 1.19 2.45 LOWER-UPSTREAM US NA NA NA
1 7/7/1999 1.59 2.47 LOWER-UPSTREAM US NA NA NA
1 7/14/1999 1.75 2.73 LOWER-UPSTREAM US NA NA NA
1 7/21/1999 1.75 2.33 LOWER-UPSTREAM US NA NA NA
1 7/28/1999 1.78 2.63 LOWER-UPSTREAM US NA NA NA
1 8/4/1999 1.37 2.25 LOWER-UPSTREAM US NA NA NA
1 8/12/1999 0.08 2.37 LOWER-UPSTREAM US 2.858553 NA 0.3560363
1 8/18/1999 2.29 2.48 LOWER-UPSTREAM US NA NA NA
1 9/4/1999 0.02 2.48 LOWER-UPSTREAM US NA NA NA
1 9/16/1999 1.93 2.45 LOWER-UPSTREAM US NA NA NA
1 10/1/1999 2.01 2.49 LOWER-UPSTREAM US NA NA NA
1 11/13/1999 1.91 2.32 LOWER-UPSTREAM US NA NA NA
1 1/4/2000 1.83 2.59 LOWER-UPSTREAM US NA NA NA
1 2/24/2000 1.93 2.83 LOWER-UPSTREAM US NA NA NA
1 3/25/2000 1.73 2.71 LOWER-UPSTREAM US NA NA NA
1 4/15/2000 1.47 3.12 LOWER-UPSTREAM US NA NA NA
1 4/28/2000 1.98 2.76 LOWER-UPSTREAM US NA NA NA
1 5/17/2000 1.80 2.92 LOWER-UPSTREAM US NA NA NA
1 5/30/2000 NA NA <NA> US 3.064330 NA 0.3970829
1 6/1/2000 1.58 2.83 LOWER-UPSTREAM US NA NA NA
1 6/6/2000 1.65 2.97 LOWER-UPSTREAM US NA NA NA
1 6/15/2000 1.85 2.65 LOWER-UPSTREAM US NA NA NA
1 6/23/2000 1.90 2.78 LOWER-UPSTREAM US NA NA NA
1 6/29/2000 1.80 2.67 LOWER-UPSTREAM US NA NA NA
1 7/6/2000 2.00 2.69 LOWER-UPSTREAM US NA NA NA
1 7/14/2000 1.97 3.46 LOWER-UPSTREAM US NA NA NA
1 7/20/2000 1.69 2.57 LOWER-UPSTREAM US NA NA NA
1 7/27/2000 1.95 3.16 LOWER-UPSTREAM US NA NA NA
1 8/2/2000 2.03 3.12 LOWER-UPSTREAM US NA NA NA
1 8/8/2000 NA NA <NA> US 3.185567 2.522467 0.3765025
1 8/11/2000 2.00 2.65 LOWER-UPSTREAM US NA NA NA
1 8/19/2000 1.95 2.50 LOWER-UPSTREAM US NA NA NA
1 9/29/2000 2.11 2.11 LOWER-UPSTREAM US NA NA NA
1 10/20/2000 2.16 2.75 LOWER-UPSTREAM US NA NA NA
1 12/1/2000 1.83 2.81 LOWER-UPSTREAM US NA NA NA
1 1/19/2001 2.01 3.01 LOWER-UPSTREAM US NA NA NA
1 3/5/2001 2.01 2.35 LOWER-UPSTREAM US NA NA NA
1 3/30/2001 2.13 2.10 LOWER-UPSTREAM US NA NA NA
1 4/21/2001 1.87 2.36 LOWER-UPSTREAM US NA NA NA
1 5/15/2001 2.00 2.37 LOWER-UPSTREAM US NA NA NA
1 5/29/2001 1.60 2.59 LOWER-UPSTREAM US NA NA NA
1 6/4/2001 NA NA <NA> US 4.233636 NA 1.3738540
1 6/7/2001 1.71 2.55 LOWER-UPSTREAM US NA NA NA
1 6/14/2001 1.83 2.75 LOWER-UPSTREAM US NA NA NA
1 6/20/2001 1.93 3.02 LOWER-UPSTREAM US NA NA NA
1 6/29/2001 1.80 2.74 LOWER-UPSTREAM US NA NA NA
1 7/6/2001 1.91 2.70 LOWER-UPSTREAM US NA NA NA
1 7/12/2001 2.04 2.46 LOWER-UPSTREAM US NA NA NA
1 7/18/2001 1.87 2.55 LOWER-UPSTREAM US NA NA NA
1 7/25/2001 1.96 2.67 LOWER-UPSTREAM US NA NA NA
1 8/2/2001 2.02 2.42 LOWER-UPSTREAM US NA NA NA
1 8/8/2001 NA NA <NA> US 2.946731 NA 0.3748240
1 8/10/2001 1.94 2.36 LOWER-UPSTREAM US NA NA NA
1 8/17/2001 1.94 0.75 LOWER-UPSTREAM US NA NA NA
1 9/1/2001 2.12 2.22 LOWER-UPSTREAM US NA NA NA
1 9/15/2001 2.15 2.35 LOWER-UPSTREAM US NA NA NA
1 9/29/2001 1.95 2.20 LOWER-UPSTREAM US NA NA NA
1 10/27/2001 1.08 1.72 LOWER-UPSTREAM US NA NA NA
1 11/17/2001 1.99 2.21 LOWER-UPSTREAM US NA NA NA
1 1/10/2002 2.01 2.31 LOWER-UPSTREAM US NA NA NA
1 2/16/2002 1.94 1.74 LOWER-UPSTREAM US NA NA NA
1 3/23/2002 2.03 2.46 LOWER-UPSTREAM US NA NA NA
1 4/13/2002 2.11 2.50 LOWER-UPSTREAM US NA NA NA
1 5/11/2002 2.11 2.54 LOWER-UPSTREAM US NA NA NA
1 5/21/2002 NA NA <NA> US 2.368000 NA 0.3235852
1 5/30/2002 2.07 2.56 LOWER-UPSTREAM US NA NA NA
1 6/6/2002 2.11 2.55 LOWER-UPSTREAM US NA NA NA
1 6/13/2002 2.03 2.66 LOWER-UPSTREAM US NA NA NA
1 6/27/2002 2.19 2.55 LOWER-UPSTREAM US NA NA NA
1 7/3/2002 2.22 2.68 LOWER-UPSTREAM US NA NA NA
1 7/11/2002 2.06 2.61 LOWER-UPSTREAM US NA NA NA
1 7/18/2002 2.23 2.58 LOWER-UPSTREAM US NA NA NA
1 7/25/2002 2.23 2.62 LOWER-UPSTREAM US NA NA NA
1 8/1/2002 2.21 2.60 LOWER-UPSTREAM US NA NA NA
1 8/8/2002 2.29 2.60 LOWER-UPSTREAM US NA NA NA
1 8/14/2002 NA NA <NA> US 1.863548 NA 0.2400621
1 8/15/2002 2.26 2.61 LOWER-UPSTREAM US NA NA NA
1 8/21/2002 2.34 2.51 LOWER-UPSTREAM US NA NA NA
1 9/6/2002 2.31 2.51 LOWER-UPSTREAM US NA NA NA
1 9/21/2002 2.25 2.55 LOWER-UPSTREAM US NA NA NA
1 10/5/2002 2.21 2.57 LOWER-UPSTREAM US NA NA NA
1 11/8/2002 2.09 2.50 LOWER-UPSTREAM US NA NA NA
1 12/23/2002 2.19 2.16 LOWER-UPSTREAM US NA NA NA
1 1/16/2003 2.23 2.15 LOWER-UPSTREAM US NA NA NA
1 3/1/2003 2.12 2.11 LOWER-UPSTREAM US NA NA NA
1 4/5/2003 2.03 1.99 LOWER-UPSTREAM US NA NA NA
1 5/2/2003 2.00 2.42 LOWER-UPSTREAM US NA NA NA
1 5/30/2003 1.66 2.09 LOWER-UPSTREAM US NA NA NA
1 6/3/2003 NA NA <NA> US 3.066333 NA 0.2926846
1 6/6/2003 1.43 2.24 LOWER-UPSTREAM US NA NA NA
1 6/12/2003 1.73 2.81 LOWER-UPSTREAM US NA NA NA
1 6/19/2003 1.51 2.97 LOWER-UPSTREAM US NA NA NA
1 7/3/2003 1.84 3.35 LOWER-UPSTREAM US NA NA NA
1 7/11/2003 2.00 2.52 LOWER-UPSTREAM US NA NA NA
1 7/17/2003 2.04 2.76 LOWER-UPSTREAM US NA NA NA
1 7/24/2003 2.02 2.40 LOWER-UPSTREAM US NA NA NA
1 7/31/2003 1.96 2.14 LOWER-UPSTREAM US NA NA NA
1 8/16/2003 2.29 1.85 LOWER-UPSTREAM US NA NA NA
1 8/20/2003 2.14 2.11 LOWER-UPSTREAM US NA NA NA
1 9/5/2003 2.03 1.45 LOWER-UPSTREAM US NA NA NA
1 10/10/2003 2.10 2.21 LOWER-UPSTREAM US NA NA NA
1 11/15/2003 1.95 1.83 LOWER-UPSTREAM US NA NA NA
1 12/22/2003 2.13 NA LOWER-UPSTREAM US NA NA NA
1 1/16/2004 2.16 2.30 LOWER-UPSTREAM US NA NA NA
1 2/13/2004 2.12 2.13 LOWER-UPSTREAM US NA NA NA
1 3/19/2004 2.02 2.11 LOWER-UPSTREAM US NA NA NA
1 4/23/2004 2.09 2.10 LOWER-UPSTREAM US NA NA NA
1 5/3/2004 1.89 2.42 LOWER-UPSTREAM US NA NA NA
1 5/17/2004 2.02 2.31 LOWER-UPSTREAM US NA NA NA
1 5/24/2004 1.97 2.58 LOWER-UPSTREAM US NA NA NA
1 5/26/2004 NA NA <NA> US 3.532614 NA 0.6329923
1 6/7/2004 1.81 2.70 LOWER-UPSTREAM US NA NA NA
1 6/14/2004 1.66 2.35 LOWER-UPSTREAM US NA NA NA
1 6/28/2004 2.10 2.70 LOWER-UPSTREAM US NA NA NA
1 7/5/2004 1.86 3.25 LOWER-UPSTREAM US NA NA NA
1 7/12/2004 2.00 2.93 LOWER-UPSTREAM US NA NA NA
1 7/19/2004 1.63 2.55 LOWER-UPSTREAM US NA NA NA
1 7/26/2004 2.01 2.65 LOWER-UPSTREAM US NA NA NA
1 8/2/2004 1.91 2.89 LOWER-UPSTREAM US NA NA NA
1 8/9/2004 2.00 2.87 LOWER-UPSTREAM US NA NA NA
1 8/17/2004 2.07 2.77 LOWER-UPSTREAM US NA NA NA
1 9/18/2004 2.18 2.67 LOWER-UPSTREAM US NA NA NA
1 10/23/2004 2.14 2.43 LOWER-UPSTREAM US NA NA NA
1 11/22/2004 2.17 2.38 LOWER-UPSTREAM US NA NA NA
1 12/20/2004 0.32 2.48 LOWER-UPSTREAM US NA NA NA
1 1/11/2005 2.16 2.44 LOWER-UPSTREAM US NA NA NA
1 2/18/2005 2.22 2.27 LOWER-UPSTREAM US NA NA NA
1 3/16/2005 2.06 2.59 LOWER-UPSTREAM US NA NA NA
1 4/9/2005 2.02 2.50 LOWER-UPSTREAM US NA NA NA
1 5/13/2005 1.96 2.62 LOWER-UPSTREAM US NA NA NA
1 6/7/2005 1.96 2.99 LOWER-UPSTREAM US NA NA NA
1 6/13/2005 1.77 2.22 LOWER-UPSTREAM US NA NA NA
1 6/20/2005 1.87 3.05 LOWER-UPSTREAM US NA NA NA
1 6/28/2005 1.76 3.50 LOWER-UPSTREAM US NA NA NA
1 7/5/2005 1.58 3.56 LOWER-UPSTREAM US NA NA NA
1 7/11/2005 1.72 3.22 LOWER-UPSTREAM US NA NA NA
1 7/19/2005 1.98 3.11 LOWER-UPSTREAM US NA NA NA
1 8/1/2005 2.14 3.02 LOWER-UPSTREAM US NA NA NA
1 8/10/2005 NA NA <NA> US 3.465650 4.032000 0.3743813
1 8/16/2005 1.79 3.00 LOWER-UPSTREAM US NA NA NA
1 10/21/2005 4.34 1.07 LOWER-UPSTREAM US NA NA NA
1 11/23/2005 2.08 1.08 LOWER-UPSTREAM US NA NA NA
1 12/21/2005 1.64 2.54 LOWER-UPSTREAM US NA NA NA
1 1/13/2006 2.05 2.01 LOWER-UPSTREAM US NA NA NA
1 2/9/2006 2.10 1.73 LOWER-UPSTREAM US NA NA NA
1 3/13/2006 2.03 2.22 LOWER-UPSTREAM US NA NA NA
1 5/22/2006 1.62 2.38 LOWER-UPSTREAM US NA NA NA
1 5/29/2006 1.61 2.29 LOWER-UPSTREAM US NA NA NA
1 6/13/2006 1.53 2.39 LOWER-UPSTREAM US NA NA NA
1 6/19/2006 1.77 2.49 LOWER-UPSTREAM US NA NA NA
1 6/20/2006 1.86 NA LOWER-UPSTREAM US NA NA NA
1 6/26/2006 1.89 3.12 LOWER-UPSTREAM US NA NA NA
1 7/3/2006 1.90 2.95 LOWER-UPSTREAM US NA NA NA
1 7/11/2006 1.26 2.49 LOWER-UPSTREAM US NA NA NA
1 7/18/2006 1.83 2.86 LOWER-UPSTREAM US NA NA NA
1 7/31/2006 1.92 2.89 LOWER-UPSTREAM US NA NA NA
1 8/15/2006 2.00 1.87 LOWER-UPSTREAM US NA NA NA
1 9/22/2006 2.07 1.62 LOWER-UPSTREAM US NA NA NA
1 11/5/2006 1.92 1.67 LOWER-UPSTREAM US NA NA NA
1 12/18/2006 2.03 1.66 LOWER-UPSTREAM US NA NA NA
1 5/24/2007 1.49 NA LOWER-UPSTREAM US NA NA NA
1 10/6/2007 2.13 NA LOWER-UPSTREAM US NA NA NA
1 12/16/2007 2.03 1.36 LOWER-UPSTREAM US NA NA NA
1 1/18/2008 2.07 1.73 LOWER-UPSTREAM US NA NA NA
1 4/19/2008 1.92 1.91 LOWER-UPSTREAM US NA NA NA
1 3/18/2010 0.61 3.20 LOWER-UPSTREAM US NA NA NA
1 5/24/2010 2.04 2.22 LOWER-UPSTREAM US NA NA NA"),header=T)
closeAllConnections()

#interpolate missing WTD values
temp.dat$WTD<-na.approx(temp.dat$WTD,as.numeric(as.Date(temp.dat$Meas_Date,"%m/%d/%Y")))

#the current limits of the left and right y-axes
ylim_1<-rev(c(0,max(na.omit(temp.dat$WTD))+1))
ylim_2<-c(floor(min(range(na.omit(temp.dat$Avg.EM.Survey.Value))[1],range(na.omit(temp.dat$Avg.Soil.Paste.EC))[1])),ceiling(max(range(na.omit(temp.dat$Avg.EM.Survey.Value))[2],range(na.omit(temp.dat$Avg.Soil.Paste.EC))[2])))

plot(as.Date(temp.dat$Meas_Date,"%m/%d/%Y"),temp.dat$WTD,xlim=c(as.Date("1999-01-01"),as.Date("2005-11-01")),ylim=ylim_1,pch=19,lwd=1.5,col="blue",cex=0.5,xlab="",ylab="",xaxt="n",type="p",las=1)
par(new=T)
plot(as.Date(temp.dat$Meas_Date,"%m/%d/%Y"),temp.dat$Avg.EM.Survey.Value,xlim=c(as.Date("1999-01-01"),as.Date("2005-11-01")),ylim=ylim_2,xlab="",ylab="WTD
(m)",xaxt="n",yaxt="n",pch=15,col="darkgreen",cex=0.75)
axis(side=4,at=seq(floor(min(range(na.omit(temp.dat$Avg.EM.Survey.Value))[1],range(na.omit(temp.dat$Avg.Soil.Paste.EC))[1])),ceiling(max(range(na.omit(temp.dat$Avg.EM.Survey.Value))[2],range(na.omit(temp.dat$Avg.Soil.Paste.EC))[2])),by=1),labels=seq(floor(min(range(na.omit(temp.dat$Avg.EM.Survey.Value))[1],range(na.omit(temp.dat$Avg.Soil.Paste.EC))[1])),ceiling(max(range(na.omit(temp.dat$Avg.EM.Survey.Value))[2],range(na.omit(temp.dat$Avg.Soil.Paste.EC))[2])),by=1),las=1)


f<-function(ylim_2){
  y2_l<-ylim_2[1]
  y2_u<-ylim_2[2]
  #solve for "x", the coeffs to map the right 
  #referenced values to their left-axis equivalent
  A<-matrix(c(y2_l,y2_u,1,1),nrow=2,ncol=2,byrow=F)
  RHS<-matrix(ylim_1,nrow=2,ncol=1,byrow=T)
  x<-solve(A,RHS)

  #map the right-referenced values to their 
  #corresponding left-axis values
  temp.dat.mapped<-x[1,1]*temp.dat$Avg.EM.Survey.Value+x[2,1]

  #calculate the resids between green squares and blue circles
 
resids<-na.omit(data.frame(Col1=temp.dat$WTD,Col2=temp.dat.mapped,SumSqrDiff=(temp.dat$WTD-temp.dat.mapped)^2))

  #quantity to be minimized
  min_this<-sum(resids$SumSqrDiff)
  min_this
}
edm<-optim(ylim_2,f,lower=c(0,0),upper=c(100,100),method="L-BFGS-B")
ylim_2<-edm$par
Tom Purucker
# 
check on:
class(temp.dat$WTD)
before and after the na.approx statement, you have changed it from
numeric to character
Best,
Tom
On Thu, Dec 30, 2010 at 3:16 PM, emorway <emorway at engr.colostate.edu> wrote:
Morway, Eric
# 
Hi Tom, 

I followed you suggestion and found the following:

class(temp.dat$WTD)
#[1] "numeric"
temp.dat$WTD<-na.approx(temp.dat$WTD,as.numeric(as.Date(temp.dat$Meas_Date,"%m/%d/%Y")))
class(temp.dat$WTD)
#[1] "numeric"

Have I misunderstood you?  If not, then my original confusion still
stands...

Respectfully,
Eric
Berend Hasselman
# 
The error occurs in the x <- solve(A,RHS) statement.
In a certain iteration ylim_2 passed by optim to your function is c(0,0)
That give a zero column of matrix A ==> exactly singular.

If I set lower=c(0,0.01) in the call of optim and use the na.approx then the
result is

ylim_2=c(0, 4.63) with the objective function=7.55

And with NO na.approx and the same lower bound the result is

ylim_2=c(0, 3.61) with the objective function=1.84

Large difference.

If you insert a print(ylim_2) in your function f you will see that after a
few iterations optim tries (0,0) for ylim_2 (when you use the na.approx).


I'm not sure if this is caused by the particular constellation of your data
or the loss of rank in your data matrix when you apply a linear
interpolation.

I have also seen that with NO na.approx optim tries ylim_2=c(0.00000000,
0.08518535) which is pretty close to 0,0).


Berend
Morway, Eric
# 
Thanks Berend, setting lower=c(0,0.01) seems to have solved the problem.  I
have a dataset with 200+ cropped fields with associated data and wanted to
generate each plot in a way that minimized the visual residuals between the
green squares and blue circles.  I have inserted a resulting plot so you can
see how your suggestion has helped me accomplish this end.

Eric

http://r.789695.n4.nabble.com/file/n3169658/Field__18_.png