How does Matlab calculate the solution of underdetermined triangular system?

I have a 4x6 Matrix R:

   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 vector b1:

b1 =

     -0.9661
      0.1590
     -0.0391
     -0.2491

To solve

R*x=b1, 
we simply do: 

x=R\b1
and the solution is:
x =

    1.0718
   -0.6078
         0
    0.2899
         0
   -0.0227

R is upper triangular, if it is square the solution is straightforward obtained by back-substitution. In the case of rectangular R, how does Matlab calculate it?

Answer 1

[Walter Roberson]

http://www.mathworks.com/help/matlab/ref/mldivide.html#zmw57dd0e481528

"Is A square?" NO -> Use QR solver

Answer 2

[John D'Errico]

That it is triangular is COMPLETELY irrelevant. In fact, I can give you a simple such system where it would be an outright terrible idea to use the triangular part directly.

For example

A = [eps*eye(4,4) + triu(rand(4,4),1),rand(4,4)]
A =
 2.2204e-16      0.43115      0.51348      0.63677      0.45466      0.65016      0.89695      0.20127
          0   2.2204e-16      0.27128      0.29868      0.46329      0.18729      0.19662      0.27423
          0            0   2.2204e-16     0.028391     0.072685       0.6784      0.94037      0.84678
          0            0            0   2.2204e-16      0.13258      0.18483      0.44909      0.97048

Here the column permuted QR solver will have no problem at all solving it, such that the problem is well conditioned. This problem is indeed well conditioned. But if we used that first triangle, we would have serious problems.

cond(A)
ans =
     10.601

If we look only at that triangular part however, it gets nasty.

cond(A(:,1:4))
ans =
   1.4785e+47

Lets try it.

b = rand(4,1);
x= A\b
x =
          0
    -1.3973
          0
          0
   0.087104
          0
    0.59965
     0.1294

As you can see, MATLAB effectively chose 4 of the unknowns to set to zero.

[Q,R,P] = qr(A,'vector');
P =
   7     8     5     2     3     6     1     4

The non-zero elements are elements [7 8 5 2], i..e, the first 4 elements of P.