Search code examples
rmathsystemlinear-algebranumerical-computing

Solving non-square linear system with R


How to solve a non-square linear system with R : A X = B ?

(in the case the system has no solution or infinitely many solutions)

Example :

A=matrix(c(0,1,-2,3,5,-3,1,-2,5,-2,-1,1),3,4,T)
B=matrix(c(-17,28,11),3,1,T)

A
     [,1] [,2] [,3] [,4]
[1,]    0    1   -2    3
[2,]    5   -3    1   -2
[3,]    5   -2   -1    1


B
     [,1]
[1,]  -17
[2,]   28
[3,]   11

Solution

  • If the matrix A has more rows than columns, then you should use least squares fit.

    If the matrix A has fewer rows than columns, then you should perform singular value decomposition. Each algorithm does the best it can to give you a solution by using assumptions.

    Here's a link that shows how to use SVD as a solver:

    http://www.ecse.rpi.edu/~qji/CV/svd_review.pdf

    Let's apply it to your problem and see if it works:

    Your input matrix A and known RHS vector B:

    > A=matrix(c(0,1,-2,3,5,-3,1,-2,5,-2,-1,1),3,4,T)
    > B=matrix(c(-17,28,11),3,1,T)
    > A
         [,1] [,2] [,3] [,4]
    [1,]    0    1   -2    3
    [2,]    5   -3    1   -2
    [3,]    5   -2   -1    1
    > B
         [,1]
    [1,]  -17
    [2,]   28
    [3,]   11
    

    Let's decompose your A matrix:

    > asvd = svd(A)
    > asvd
    $d
    [1] 8.007081e+00 4.459446e+00 4.022656e-16
    
    $u
               [,1]       [,2]       [,3]
    [1,] -0.1295469 -0.8061540  0.5773503
    [2,]  0.7629233  0.2908861  0.5773503
    [3,]  0.6333764 -0.5152679 -0.5773503
    
    $v
                [,1]       [,2]       [,3]
    [1,]  0.87191556 -0.2515803 -0.1764323
    [2,] -0.46022634 -0.1453716 -0.4694190
    [3,]  0.04853711  0.5423235  0.6394484
    [4,] -0.15999723 -0.7883272  0.5827720
    
    > adiag = diag(1/asvd$d)
    > adiag
              [,1]      [,2]        [,3]
    [1,] 0.1248895 0.0000000 0.00000e+00
    [2,] 0.0000000 0.2242431 0.00000e+00
    [3,] 0.0000000 0.0000000 2.48592e+15
    

    Here's the key: the third eigenvalue in d is very small; conversely, the diagonal element in adiag is very large. Before solving, set it equal to zero:

    > adiag[3,3] = 0
    > adiag
              [,1]      [,2] [,3]
    [1,] 0.1248895 0.0000000    0
    [2,] 0.0000000 0.2242431    0
    [3,] 0.0000000 0.0000000    0
    

    Now let's compute the solution (see slide 16 in the link I gave you above):

    > solution = asvd$v %*% adiag %*% t(asvd$u) %*% B
    > solution
              [,1]
    [1,]  2.411765
    [2,] -2.282353
    [3,]  2.152941
    [4,] -3.470588
    

    Now that we have a solution, let's substitute it back to see if it gives us the same B:

    > check = A %*% solution
    > check
         [,1]
    [1,]  -17
    [2,]   28
    [3,]   11
    

    That's the B side you started with, so I think we're good.

    Here's another nice SVD discussion from AMS:

    http://www.ams.org/samplings/feature-column/fcarc-svd