CODEWITHSUNDEEP
rmatlab
Sinem
Sinem

Reputation: 63

How does Matlab calculate the solution of underdetermined systems?

For solving rectangular linear systems Ax=b, where A is mxn and n>m, Matlab performs first a QR factorization then solves a triangular linear system by substitution. For example: if A is a 4x6 matrix:

A =

    0.8147    0.9058    0.1270    0.9134    0.6324    0.0975
    0.2785    0.5469    0.9575    0.9649    0.1576    0.9706
    0.9572    0.4854    0.8003    0.1419    0.4218    0.9157
    0.7922    0.9595    0.6557    0.0357    0.8491    0.9340

and 

b=

   -0.9661
    0.1590
   -0.0391
   -0.2491

to solve Ax=b we simply do:

x=A\b;

and the solution is:

x =

   -0.4284
   -0.6475
         0
   -0.1153
         0
    0.7662

The underlying steps in this computation are:

  1. QR decomposition of A, then Ax=bwould be equivalent to Q*R*x=b.
  2. As inv(Q)=Q'and R is upper triangular, solving the system would return to solving R*x=b1 where b1=Q'*b.

As R is upper triangular of size 6x4, we have to perform a back-substitution. How does Matlab perform the back substitution knowing that R is not square?

I'm porting this computation to C and have succesfully done it until R*x=b1, I'm lost with the back substitution and I want to find the same result as Matlab.

Edit: For the above matrix A: Solving Ax=bis equivalent to solving Rx=b1where:

R =

   -1.5117   -1.3991   -1.0952   -0.7786   -1.0819   -1.3007
         0   -0.5641   -0.2197   -0.6538   -0.2920   -0.2481
         0         0   -0.8692   -0.2077    0.1422   -0.9295
         0         0         0   -0.8426    0.2182    0.2125

and 

b1 =

   -0.9661
    0.1590
   -0.0391
   -0.2491

This system is solved by back-substitution as R is upper triangular, this is straightforward in case of square triangular matrices (i.e m=n), However, R is not square. How does Matlab perform the back-substitution to find x?

Upvotes: 2

Views: 1385

Answers (2)

pragmatist1
pragmatist1

Reputation: 919

So I'm not 100% sure what MATLAB does, but I can tell you the following:

For your underdetermined matrix A (nxm), we can take A' = QR. We also see that A = R'Q'. Note that now, R has n nonzero rows. We can then solve the system by recognizing R'Q'x=b, then Q'x=inv(R')*b. The RHS can be solved by back-substitution, we'll call it y. So Q'x=y. But Q is self-Hermitian, so x=Qy. And there you have the least squares solution.

Doing this in MATLAB we get:

>> [Q,R] = qr(A')

Q =

   -0.4918    0.2143   -0.6131   -0.5675   -0.1086    0.0503
   -0.5468    0.0638    0.0596    0.5238   -0.5922   -0.2614
   -0.0767   -0.6389   -0.1919   -0.0733    0.2275   -0.7014
   -0.5514   -0.2397    0.6743   -0.3550    0.1568    0.1821
   -0.3818    0.2094   -0.1834    0.4750    0.7368    0.0904
   -0.0589   -0.6637   -0.3090    0.2158   -0.1351    0.6291


R =

   -1.6565   -1.1588   -1.0907   -1.3634
         0   -1.3597   -0.8286   -0.6385
         0         0   -0.9760   -0.9745
         0         0         0    0.5972
         0         0         0         0
         0         0         0         0

>> Q*(R'\b)

ans =

   -0.4256
   -0.3057
    0.3568
   -0.2745
   -0.2823
    0.4249

>>

I used backslash above for R'\b, but it's clear you can back-substitute in your own routine.

We can verify this by computing the solution using the pseudo-inverse of A:

>> A'*inv(A*A')*b

ans =

   -0.4256
   -0.3057
    0.3568
   -0.2745
   -0.2823
    0.4249

You can easily check that this is a valid solution of course. This should be readily implementable in C.

Hope this helps.

Upvotes: 2

gregswiss
gregswiss

Reputation: 1456

How MATLAB implements in details the under-determined case only Mathworks know! But a way you can do it is as follow.

A*x = b
A' = Q*R

Now if you look at R you will have something of the form

R = [R1; zeros]

where R1 is square.

you can express the solution as:

x = Q * [inv(R1')*b; zeros]

as you indicates that you already solved the problem of inverting a square matrix that should be easy :)

Upvotes: 1

Related Questions