Here's a simple example to illustrate what I want:

I have the following linear system:

2*x + 1*y + 0*z + 0*w = 3
1*x  - 1*y + 0*z + 0*w = 0
0*x + 0*y + 1*z + 1*w = 1
1*x - 1*y + 2*z + 2*w = 2

which can be written in the form:  A.x = b

The system is underdetermined:  rank (A) = 3

I can solve for "x" and "y" but "z" and "w" are linearly dependent
variables.
I can only find a complete solution by setting a value for either "z"
or "w"
(or add another independent relationship between the two, but I want
to work with the "original" system. I also dont want to put additional
constraint like minimal norm solutions)

So here's my question:

For any linear system ("A" not necessarily square), how can I know
which variables are readily calculable and then get a suggestion for
the variables that i need to specify for solving the system
(suggestion in the sense that there is not a unique answer).

This is part of an optimization routine where I need to find values of
"x", "y", "z", "w" .... satisfying a certain objective function, but I
dont the optimization routine to:
1) Not change the invariant ("x" and "y")
2) Optimize only the minimum number of unknown ("w" or "z", but not
both).

hope it's clear.

little typo:

> This is part of an optimization routine where I need to find values of
> "x", "y", "z", "w" .... satisfying a certain objective function, but I
> WANT the optimization routine to:
> 1) Not change the invariant ("x" and "y")
> 2) Optimize only the minimum number of unknown ("w" or "z", but not
> both).
>
> hope it's clear.

On Jan 6, 4:21=A0pm, Euh <Euh...@hotmail.com> wrote:
> Here's a simple example to illustrate what I want:
>
> I have the following linear system:
>
> 2*x + 1*y + 0*z + 0*w =3D 3
> 1*x =A0- 1*y + 0*z + 0*w =3D 0
> 0*x + 0*y + 1*z + 1*w =3D 1
> 1*x - 1*y + 2*z + 2*w =3D 2
>
> which can be written in the form: =A0A.x =3D b
>
> The system is underdetermined: =A0rank (A) =3D 3
>
> I can solve for "x" and "y" but "z" and "w" are linearly dependent
> variables.
> I can only find a complete solution by setting a value for either "z"
> or "w"
> (or add another independent relationship between the two, but I want
> to work with the "original" system. I also dont want to put additional
> constraint like minimal norm solutions)
>
> So here's my question:
>
> For any linear system ("A" not necessarily square), how can I know
> which variables are readily calculable and then get a suggestion for
> the variables that i need to specify for solving the system
> (suggestion in the sense that there is not a unique answer).
>
> This is part of an optimization routine where I need to find values of
> "x", "y", "z", "w" .... satisfying a certain objective function, but I
> dont the optimization routine to:
> 1) Not change the invariant ("x" and "y")
> 2) Optimize only the minimum number of unknown ("w" or "z", but not
> both).
>
> hope it's clear.

nullA =3D
         0
   -0.0000
   -0.7071
    0.7071

Therefore each solution has the same x, y, and w+z.

It is obvious if you stare at your equations long enough.

Hope this helps.

Greg

On Jan 6, 7:53=A0pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 6, 4:21=A0pm, Euh <Euh...@hotmail.com> wrote:
>
>
>
>
>
> > Here's a simple example to illustrate what I want:
>
> > I have the following linear system:
>
> > 2*x + 1*y + 0*z + 0*w =3D 3
> > 1*x =A0- 1*y + 0*z + 0*w =3D 0
> > 0*x + 0*y + 1*z + 1*w =3D 1
> > 1*x - 1*y + 2*z + 2*w =3D 2
>
> > which can be written in the form: =A0A.x =3D b
>
> > The system is underdetermined: =A0rank (A) =3D 3
>
> > I can solve for "x" and "y" but "z" and "w" are linearly dependent
> > variables.
> > I can only find a complete solution by setting a value for either "z"
> > or "w"
> > (or add another independent relationship between the two, but I want
> > to work with the "original" system. I also dont want to put additional
> > constraint like minimal norm solutions)
>
> > So here's my question:
>
> > For any linear system ("A" not necessarily square), how can I know
> > which variables are readily calculable and then get a suggestion for
> > the variables that i need to specify for solving the system
> > (suggestion in the sense that there is not a unique answer).
>
> > This is part of an optimization routine where I need to find values of
> > "x", "y", "z", "w" .... satisfying a certain objective function, but I
> > dont the optimization routine to:
> > 1) Not change the invariant ("x" and "y")
> > 2) Optimize only the minimum number of unknown ("w" or "z", but not
> > both).
>
> > hope it's clear.
>
> nullA =3D
> =A0 =A0 =A0 =A0 =A00
> =A0 =A0-0.0000
> =A0 =A0-0.7071
> =A0 =A0 0.7071
>
> Therefore each solution has the same x, y, and w+z.
>
> It is obvious if you stare at your equations long enough.

Since the system is undetermined,

x =3D pinv(A)*b   is a solution (minimum norm).

 t
To get another solution add any constant times the null solution.


Hope this helps.

Greg-

On 6 jan, 19:53, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 6, 4:21=A0pm, Euh <Euh...@hotmail.com> wrote:
>
>
>
> > Here's a simple example to illustrate what I want:
>
> > I have the following linear system:
>
> > 2*x + 1*y + 0*z + 0*w =3D 3
> > 1*x =A0- 1*y + 0*z + 0*w =3D 0
> > 0*x + 0*y + 1*z + 1*w =3D 1
> > 1*x - 1*y + 2*z + 2*w =3D 2
>
> > which can be written in the form: =A0A.x =3D b
>
> > The system is underdetermined: =A0rank (A) =3D 3
>
> > I can solve for "x" and "y" but "z" and "w" are linearly dependent
> > variables.
> > I can only find a complete solution by setting a value for either "z"
> > or "w"
> > (or add another independent relationship between the two, but I want
> > to work with the "original" system. I also dont want to put additional
> > constraint like minimal norm solutions)
>
> > So here's my question:
>
> > For any linear system ("A" not necessarily square), how can I know
> > which variables are readily calculable and then get a suggestion for
> > the variables that i need to specify for solving the system
> > (suggestion in the sense that there is not a unique answer).
>
> > This is part of an optimization routine where I need to find values of
> > "x", "y", "z", "w" .... satisfying a certain objective function, but I
> > dont the optimization routine to:
> > 1) Not change the invariant ("x" and "y")
> > 2) Optimize only the minimum number of unknown ("w" or "z", but not
> > both).
>
> > hope it's clear.
>
> nullA =3D
> =A0 =A0 =A0 =A0 =A00
> =A0 =A0-0.0000
> =A0 =A0-0.7071
> =A0 =A0 0.7071
>
> Therefore each solution has the same x, y, and w+z.
>
> It is obvious if you stare at your equations long enough.
>

Yup that was a simple system for illustration. But with a bigger
matrix, how do i interpret the null space ?

- lines with zeros are calculable ?
- how do you "pair" the remaining dependent variables ? (in my example
it was obvious, but for the general case ?)

Euh <EuhSCQ@hotmail.com> wrote in message <88fd67ce-82d9-40a1-bddd-67e88a21865b@o28g2000yqh.googlegroups.com>...
>
> 
> So here's my question:
> 
> For any linear system ("A" not necessarily square), how can I know
> which variables are readily calculable and then get a suggestion for
> the variables that i need to specify for solving the system
> (suggestion in the sense that there is not a unique answer).
> 

Use QR decomposition with row-permutation of the matrix (three output form). This will tell you which non-redundant parameters you can keep to have the same rank as the original one.

This is basically along the line the slash ("\") operator does.

Bruno

Sorry typo correction:

> 
> Use QR decomposition with *column-permutation* of the matrix (three output form). This will tell you which non-redundant parameters you can keep to have the same rank as the original one.

Bruno

Euh <EuhSCQ@hotmail.com> wrote in message <88fd67ce-82d9-40a1-bddd-67e88a21865b@o28g2000yqh.googlegroups.com>...

> So here's my question:
> 
> For any linear system ("A" not necessarily square), how can I know
> which variables are readily calculable and then get a suggestion for
> the variables that i need to specify for solving the system
> (suggestion in the sense that there is not a unique answer).


Here's yet another solution using rref().

%data
A=[2 1 0 0 ; 1 -1 0 0; 0 0 1 1; 1 -1 2 2]; 
b=[3 0 1 2].';
 
%engine
[R,NonRedundantParams]=rref(A); 
Solution=zeros(size(A,2),1); 
xb=A(:,NonRedundantParams)\b;

Solution(NonRedundantParams)=xb

Solution =

    1.0000
    1.0000
    1.0000
         0

This procedure presumes, by the way, that the system has a solution at all. You haven't said how you would like to deal with the case where the solution space is empty. However, you can of course verify success by computing

>> Error=norm( A*Solution-b )

Error =

  1.1957e-015

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hi4d41$7a7$1@fred.mathworks.com>...

> 
> 
> Here's yet another solution using rref().
> 

Careful, RREF is an academic code. It is not recommended to use it beside that boundary.

Bruno 

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hi4dqd$4im$1@fred.mathworks.com>...

> 
> Careful, RREF is an academic code. It is not recommended to use it beside that boundary.

Do you mean that RREF uses an unreliable algorithm, or that it is simply slower because it is in M-code, rather than a built-in function like QR?

 

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hi4fem$h4b$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hi4dqd$4im$1@fred.mathworks.com>...
> 
> > 
> > Careful, RREF is an academic code. It is not recommended to use it beside that boundary.
> 
> Do you mean that RREF uses an unreliable algorithm, or that it is simply slower because it is in M-code, rather than a built-in function like QR?

Both.

RREF is in principle less costly in number of calculation than QR, but the rank detected is less reliable, more affected by round off error. To detect the rank, the order of method is (more reliable first): SVD then QR then RREF.

RREF is slow because it is implemented as *.m, and if you take a look at the code, there is a a "strange" fractional rounding at the end. This code is programmed by Cleve Moler between two lunches.

Bruno

On 7 jan, 03:35, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> Sorry typo correction:
>
>
>
> > Use QR decomposition with *column-permutation* of the matrix (three output form). This will tell you which non-redundant parameters you can keep to have the same rank as the original one.
>
> Bruno


let use this example:

Ax=0

with

A=[     0     0     0     0     0     0     0     0     0     1
-1     0     0
     0     0     0     0     0     0     0    -1     1     0     0
0     0
     0     1     0     0     0     0     0     1     1    -1     0
0     0
     0     0     0     0     0     0     1     1    -1     1     1
-1     0
     1    -1    -1     0     0     0     0     0     0     0     0
0     0
     0     0     0     0     0     0     0     0     0     0     0
1    -1
     0     0     0     0     1    -1    -1     0     0     0     0
0     0
     0     0     1    -1    -1     0     0     0     0     0     0
0     0
     0     0     0     0     0     0     1     0    -1     0     0
0     0
     0     0     0     1     0     0    -1    -1     0     0     0
0     0];

Assuming x(1) and x(6) are known:   x(1)=0.134 and x(6)=0


I know that by specifying just x(2), I can generate a solution:

For x(2) =0.08 I get:

 x=[0.1340
    0.0800
    0.0540
    0.0360
    0.0180
         0
    0.0180
    0.0180
    0.0180
    0.1160
    0.1160
    0.2500
    0.2500];

I'm just not sure how to generate a code that could tell me that x(2)
is sufficient given that x(1) and x(6) are known

The null space of A has 3 column. By fixing 2 variables, I know I only
have one more to specify (right ?). But which one ?

% Do this:

A=[ 0 0 0 0 0 0 0 0 0 1 -1 0 0
     0 0 0 0 0 0 0 -1 1 0 0 0 0
     0 1 0 0 0 0 0 1 1 -1 0 0 0
     0 0 0 0 0 0 1 1 -1 1 1 -1 0
     1 -1 -1 0 0 0 0 0 0 0 0 0 0
     0 0 0 0 0 0 0 0 0 0 0 1 -1
     0 0 0 0 1 -1 -1 0 0 0 0 0 0
     0 0 1 -1 -1 0 0 0 0 0 0 0 0
     0 0 0 0 0 0 1 0 -1 0 0 0 0
     0 0 0 1 0 0 -1 -1 0 0 0 0 0];

allidx = 1:size(A,2);
knownidx = [1 6];
remainedidx = setdiff(allidx,knownidx);

% Remove column 1 & 6
B = A(:,remainedidx);
[Q R P] = qr(B,0);

>> diag(R)

ans =

    2.0000
    1.7321
    1.7321
    1.4142
    1.3844
    1.2218
    1.0358
    1.0334
    0.9829
    0.1987

>> size(diag(R))

ans =

    10    1

% That means there are 10 independent columns of B: P(1:10);
% The 11th column of B = B(:,P(11))) is dependent of the rest.
%
% If diag(R) has elements "equal" to zero you should take all those indexes as well
% as dependent columns. Note: "equal" must be check with a tolerance

>> P(11) % column of B that depends on the rest

ans =

     1

>> remainedidx(P(11)) % this is the corresponding column of A

ans =

     2

% Bruno

On Jan 7, 5:30=A0am, "Matt J " <mattjacREM...@THISieee.spam> wrote:
> Euh <Euh...@hotmail.com> wrote in message <88fd67ce-82d9-40a1-bddd-67e88a=
218...@o28g2000yqh.googlegroups.com>...
> > So here's my question:
>
> > For any linear system ("A" not necessarily square), how can I know
> > which variables are readily calculable and then get a suggestion for
> > the variables that i need to specify for solving the system
> > (suggestion in the sense that there is not a unique answer).
>
> Here's yet another solution using rref().
>
> %data
> A=3D[2 1 0 0 ; 1 -1 0 0; 0 0 1 1; 1 -1 2 2];
> b=3D[3 0 1 2].';
>
> %engine
> [R,NonRedundantParams]=3Drref(A);
> Solution=3Dzeros(size(A,2),1);
> xb=3DA(:,NonRedundantParams)\b;
>
> Solution(NonRedundantParams)=3Dxb
>
> Solution =3D
>
> =A0 =A0 1.0000
> =A0 =A0 1.0000
> =A0 =A0 1.0000
> =A0 =A0 =A0 =A0 =A00
>
> This procedure presumes, by the way, that the system has a solution at al=
l. You haven't said how you would like to deal with the case where the solu=
tion space is empty.

If the system is underdetermined won't there
always be a solution?

Hope this helps.

Greg

On Jan 7, 3:35 am, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:

> Use QR decomposition with *column-permutation*
> of the matrix (three output form). This will tell
> you which non-redundant parameters you can
> keep to have the same rank as the original one.

Consider

A*x     = b, [A b} is given, x is unknown
size(A) = [m n], m <= n (Not overdetermined)
size(x) = [n 1]
size(b) = [m 1]
A'*b   ~= 0     (A solution exists)
rank(A) = r < m  (Infinite number of solutions exist)

Then the "basic" solutions contain at least
n-r zero elements

[Q R PT] = qr(A)
A = Q*R*P,   P = PT'
size(Q) = [m m]
size(R) = [m n] upper triangular
size(P) = [n n]

The last n-r rows of P indicate components
of x that can be arbitrarily specified.
For a "basic" solution they are chosen to be zero.
The first r rows of P indicate nonzero
components of a "basic" solution

For example,

close all, clear all, clc

A = [ 2  1 0 0 ;
      1 -1 0 0 ;
      0  0 1 1;
      1 -1 2 2]

 b = [3 0 1 2]'
 [n q]  = size(b)                     % [ 4 1]
 [m n]  = size(A)                   % [ 4 4]
 r      = rank(A)                     % 3
 condA  = cond(A)                % Inf
 xexist = norm(A'*b)~=0         % 1

[Q R PT] = qr(A)

[m m]    = size(Q)
[m n]    = size(R)
P        = PT'
[n n]    = size(P)
%  P =
%      1     0     0     0
%      0     0     1     0
%      0     1     0     0
%      0     0     0     1


basic    = P(1:r,:)*[1:n]'         % [ 1 3 2 ]'
nonbasic = P(r+1:n,:)*[1:n]'   % 4

Hope this helps.

Greg

On Jan 7, 10:11=A0am, Euh <Euh...@hotmail.com> wrote:
> On 7 jan, 03:35, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
>
> > Sorry typo correction:
>
> > > Use QR decomposition with *column-permutation*
> > >of the matrix (three output form). This will tell you which
> > >non-redundant parameters you can keep to have the
> > >same rank as the original one.

>
> let use this example:
>
> Ax=3D0
>
> with
>
> A=3D[ =A0...
>      0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 1  -1   0   0
>  =A0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0-1 =A0 1 =A0 0 =A0 0   =
0 =A0 0
> =A0 =A0 =A00 =A0 1 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 1 =A0 1 =A0-1 =A0 0 =
  0 =A0 0
> =A0 =A0 =A00 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 1 =A0 1 =A0-1 =A0 1 =A0 1 =
 -1 =A0 0
> =A0 =A0 =A01 =A0-1 =A0-1 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =
  0 =A0 0
> =A0 =A0 =A00 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =
  1 =A0-1
> =A0 =A0 =A00 =A0 0 =A0 0 =A0 0 =A0 1 =A0-1 =A0-1 =A0 0 =A0 0 =A0 0 =A0 0 =
  0 =A0 0
> =A0 =A0 =A00 =A0 0 =A0 1 =A0-1 =A0-1 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =
  0 =A0 0
> =A0 =A0 =A00 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 1 =A0 0 =A0-1 =A0 0 =A0 0 =
  0 =A0 0
> =A0 =A0 =A00 =A0 0 =A0 0 =A0 1 =A0 0 =A0 0 =A0-1 =A0-1 =A0 0 =A0 0 =A0 0 =
  0 =A0 0];
>
> Assuming x(1) and x(6) are known: =A0 x(1)=3D0.134 and x(6)=3D0
>
> I know that by specifying just x(2), I can generate a solution:

I don't understand your logic.

Consider this:

size(A) =3D [10 13]
rank(A) =3D 10
Therefore 13 - 10 =3D 3 components can be
set to zero for a "basic" solution.

Although columns 1, 6 and 13 have only one component,
using the QRP decomposition method of my previous post,
reveals that [x1 x6 x13] =3D 0 constitutes a "basic" solution

 [m n]  =3D size(A)               % [ 10 13]
 r      =3D rank(A)               % 10
 condA  =3D cond(A)               % 6.9666

[Q R PT] =3D qr(A)
 P       =3D PT'
[m m]    =3D size(Q)
[m n]    =3D size(R)
[n n]    =3D size(P)

P =3D [...
% 0     0     0     0     0     0     1     0     0     0     0
0     0;
% 0     0     0     0     0     0     0     1     0     0     0
0     0;
% 0     0     0     0     0     0     0     0     1     0     0
0     0;
% 0     0     1     0     0     0     0     0     0     0     0
0     0;
% 0     0     0     0     0     0     0     0     0     1     0
0     0;
% 0     0     0     0     0     0     0     0     0     0     1
0     0;
% 0     0     0     0     1     0     0     0     0     0     0
0     0;
% 0     0     0     0     0     0     0     0     0     0     0
1     0;
% 0     0     0     1     0     0     0     0     0     0     0
0     0;
% 1     0     0     0     0     0     0     0     0     0     0
0     0;
% 0     0     0     0     0     1     0     0     0     0     0
0     0;
% 0     1     0     0     0     0     0     0     0     0     0
0     0;
% 0     0     0     0     0     0     0     0     0     0     0
0     1]


basic    =3D P(1:r,:)*[1:n]'     % [ 7 8 9 3 10 11 5 12 4 1]'
nonbasic =3D P(r+1:n,:)*[1:n]'   % [6 2 13]'

Hope this helps.

Greg

On Jan 9, 1:59=A0am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 7, 10:11=A0am, Euh <Euh...@hotmail.com> wrote:
>
>
>
>
>
> > On 7 jan, 03:35, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
>
> > > Sorry typo correction:
>
> > > > Use QR decomposition with *column-permutation*
> > > >of the matrix (three output form). This will tell you which
> > > >non-redundant parameters you can keep to have the
> > > >same rank as the original one.
>
> > let use this example:
>
> > Ax=3D0
>
> > with
>
> > A=3D[ =A0...
> > =A0 =A0 =A00 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 1 =A0-=
1 =A0 0 =A0 0
> > =A0=A0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0-1 =A0 1 =A0 0 =A0 =
0 =A0 0 =A0 0
> > =A0 =A0 =A00 =A0 1 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 1 =A0 1 =A0-1 =A0 =
0 =A0 0 =A0 0
> > =A0 =A0 =A00 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 1 =A0 1 =A0-1 =A0 1 =A0 =
1 =A0-1 =A0 0
> > =A0 =A0 =A01 =A0-1 =A0-1 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 =
0 =A0 0 =A0 0
> > =A0 =A0 =A00 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 =
0 =A0 1 =A0-1
> > =A0 =A0 =A00 =A0 0 =A0 0 =A0 0 =A0 1 =A0-1 =A0-1 =A0 0 =A0 0 =A0 0 =A0 =
0 =A0 0 =A0 0
> > =A0 =A0 =A00 =A0 0 =A0 1 =A0-1 =A0-1 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 =
0 =A0 0 =A0 0
> > =A0 =A0 =A00 =A0 0 =A0 0 =A0 0 =A0 0 =A0 0 =A0 1 =A0 0 =A0-1 =A0 0 =A0 =
0 =A0 0 =A0 0
> > =A0 =A0 =A00 =A0 0 =A0 0 =A0 1 =A0 0 =A0 0 =A0-1 =A0-1 =A0 0 =A0 0 =A0 =
0 =A0 0 =A0 0];
>
> > Assuming x(1) and x(6) are known: =A0 x(1)=3D0.134 and x(6)=3D0
>
> > I know that by specifying just x(2), I can generate a solution:
>
> I don't understand your logic.
>
> Consider this:
>
> size(A) =3D [10 13]
> rank(A) =3D 10
> Therefore 13 - 10 =3D 3 components can be
> set to zero for a "basic" solution.
>
> Although columns 1, 6 and 13 have only one component,
> using the QRP decomposition method of my previous post,
> reveals that [x1 x6 x13] =3D 0 constitutes a "basic" solution

WHOOPS! That should be [x2 x6  x13] =3D 0
>
> =A0[m n] =A0=3D size(A) =A0 =A0 =A0 =A0 =A0 =A0 =A0 % [ 10 13]
> =A0r =A0 =A0 =A0=3D rank(A) =A0 =A0 =A0 =A0 =A0 =A0 =A0 % 10
> =A0condA =A0=3D cond(A) =A0 =A0 =A0 =A0 =A0 =A0 =A0 % 6.9666
>
> [Q R PT] =3D qr(A)
> =A0P =A0 =A0 =A0 =3D PT'
> [m m] =A0 =A0=3D size(Q)
> [m n] =A0 =A0=3D size(R)
> [n n] =A0 =A0=3D size(P)
>
> P =3D [...
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 1 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 1=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 1 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 1 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 1 =A0 =A0 0
> 0 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 1
> 0 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 1 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 1 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 1 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 0;
> % 1 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 1 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 0;
> % 0 =A0 =A0 1 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 0;
> % 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0=
 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
> 0 =A0 =A0 1]
>
> basic =A0 =A0=3D P(1:r,:)*[1:n]' =A0 =A0 % [ 7 8 9 3 10 11 5 12 4 1]'
> nonbasic =3D P(r+1:n,:)*[1:n]' =A0 % [6 2 13]'

 Hope this helps.

 Greg-