0

We are trying to do nonlinear least squares fitting in R. We have reference code in MATLAB as below.

 A =  [ -4.09549023 -1.5924967 -5.2775267 0.7195365 -1.4681932;
 -0.09302291  0.2538085  0.6080219 0.1413484 -0.4614447;
    1.00000000  1.0000000  1.0000000 1.0000000  1.0000000;]

 b = [-0.52701539;0.08974224;1.00000000]

[w,res] = lsqnonneg(A, b)

w =

         0
         0
    0.1377
    0.6700
    0.1922


res =

   1.2519e-32

If we try the same thing in R, the results are entirely different.

A <- rbind(cbind(4.09549023, -1.5924967, -5.2775267, 0.7195365, -1.4681932),
cbind(-0.09302291,  0.2538085,  0.6080219, 0.1413484, -0.4614447),
 cbind(1.00000000,  1.0000000,  1.0000000, 1.0000000,  1.0000000))

b <- rbind(-0.52701539,0.08974224,1.00000000)

sol <- lsqnonneg(A,as.vector(b))

$x
[1] 0.2212015 0.1986838 0.1890215 0.2080815 0.1830117

$resnorm
[1] 4.889664e-10

$residual
              [,1]
[1,]  2.871547e-11
[2,]  4.866847e-10
[3,] -3.743827e-11

The problem is, though the input values of A and b are the same, and the methods invoked both in R and MATLAB have the same purpose (viz., non-negative least squares fitting) the residuals are different in both cases. resnorm = 4.889664e-10 in R and res = 1.2519e-32 in MATLAB. If anyone can point out some obvious mistake in the results, or suggest some alternative method in R to match the method in MATLAB, it would be really appreciated.

Nick Stauner
  • 11,558
  • 5
  • 47
  • 105
Ravindran_Quan
  • 1
  • 2
  • 1
    Although I formatted this question for readability, it is still incomprehensible. Neither of the code sections looks like valid code or output--pieces evidently have been stripped out of them. How are we supposed to be comparing the results? You need to explain how you are interpreting this output and point out specific discrepancies. When you are making those edits, you might also find the varying answers and comments at http://stats.stackexchange.com/questions/7308 worth reading, because they appear to deal with a similar issue. – whuber May 26 '14 at 15:16
  • 1
    (1) If you use an R package, at least tell us which one you used. Is this the `lsqnonneg` function in `pracma` or something else? (2) Even if your solution were determined by your equations, you still wouldn't expect identical sums of squared residuals, since the two won't by default use the same stopping criteria, nor likely even the same algorithms. – Glen_b May 26 '14 at 23:59

1 Answers1

3

Your system of equations $Ax=b$ has more variables than equations and happens to be consistent. Thus there are infinitely many solutions. It appears that there are still infinitely many solutions when you include the nonnegativity constraint. This behavior isn't at all surprising- why did you expect to get the same answer from both packages when the underlying problem has infinitely many solutions? Is there some particular answer that you'd like to get, and why?

Brian Borchers
  • 5,015
  • 1
  • 18
  • 27
  • It's clear that the tolerance used by the R routine is fairly loose in comparison with MATLAB, but in both cases the norm of the residuals is essentially 0. – Brian Borchers May 27 '14 at 15:06