0

I want to find the least square estimate for an over-determined system of linear equations the way Matlab does with \.

What I am trying to reproduce in R:

% matlab code

X = [3.8642    9.6604;
    14.2000   35.5000;
    41.7832  104.4580;
     0.4084    1.0210];

y = [1.2300
     4.5200
    13.3000
     0.1300];

X\y  %  => [0, 0.1273]

I tried R's lsfit method, the generalized Inverse (ginv) from the MASS package, and using the QR compositon (R-1Q'y), but all returns different results.

Data in R format:

x <- matrix(c(3.8642, 14.2, 41.7832, 0.4084, 9.6604, 35.5, 104.458, 1.021),
            ncol = 2)
y <- c(1.23, 4.52, 13.3, 0.13)
Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
M. F. Bostanci
  • 75
  • 2
  • 5
  • 3
    Good job providing the MATLAB code you want to emulate, but in terms of a [mcve], it would be nice to also provide the data in a way for people to easily cut & paste into R – duckmayr Dec 16 '19 at 13:32

2 Answers2

3

How to do the equivalent of MATLAB's \ for least squares estimation? The documentation for \ says

x = A\B solves the system of linear equations A*x = B.

The equivalent in R you're looking for is solve() or qr.solve(), in this case:

?solve

This generic function solves the equation a %*% x = b for x
. . .
qr.solve can handle non-square systems.

So, we can see

qr.solve(x, y)
# [1] 0.003661243 0.125859408

This is very close to your MATLAB solution. Likewise, lsfit() or lm() (of course) give you the same answer:

coef(lm(y ~ x + 0))
#          x1          x2 
# 0.003661243 0.125859408 
lsfit(x, y, intercept = FALSE)$coef
#          X1          X2 
# 0.003661243 0.125859408

We can see this answer fit the data at least as well as your MATLAB solution:

r_solution <- coef(lm(y ~ x + 0))
y - (x %*% r_solution)
              [,1]
[1,] -1.110223e-15
[2,]  1.366296e-06
[3,] -4.867456e-07
[4,]  2.292817e-06
matlab_solution <- c(0, 0.1273)
y - (x %*% matlab_solution)
           [,1]
[1,] 0.00023108
[2,] 0.00085000
[3,] 0.00249660
[4,] 0.00002670
Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
duckmayr
  • 16,303
  • 3
  • 35
  • 53
1

When I used ginv(), I am able to achieve the same result as MATLAB:

  • using ginv() in R
library(MASS)
X <- matrix(c(3.8642, 14.2, 41.7832, 0.4084, 9.6604, 35.5, 104.458, 1.021),ncol = 2)
y <- matrix(c(1.23,4.52,13.3,0.13),ncol = 1)

> ginv(X)%*%y
            [,1]
[1,] 0.003661243
[2,] 0.125859408
  • using pinv() or \ in MATLAB
X = [3.8642    9.6604;
   14.2000   35.5000;
   41.7832  104.4580;
    0.4084    1.0210];

y = [1.2300
    4.5200
   13.3000
    0.1300];

>> X\y
ans =

   0.0036612
   0.1258594

>> pinv(X)*y
ans =

   0.0036612
   0.1258594
ThomasIsCoding
  • 96,636
  • 9
  • 24
  • 81
  • 1
    MATLAB is not Octave, and might not use the same algorithms under all circumstances to solve systems of linear equations. Don't assume that because you got some answer in Octave, that MATLAB would produce the same answer. – Cris Luengo Dec 07 '22 at 21:24