CODEWITHSUNDEEP
rmatlablinear-algebralinear-equation
M. F. Bostanci
M. F. Bostanci

Reputation: 87

Matlab's least square estimate via \ for over-determined system in R

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)

Upvotes: 0

Views: 256

Answers (2)

ThomasIsCoding
ThomasIsCoding

Reputation: 101628

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

Upvotes: 1

duckmayr
duckmayr

Reputation: 16930

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

Upvotes: 3

Related Questions