Basis for null space of a matrix

[Dominique Vercammen]

Hello all,

I was trying to find an orthonormal basis for the null space of an SXMatrix using QR decomposition in CasADi, but this method only returns the left part of the Q matrix, not the right part of Q which is the orthonormal basis.

Is there a way to get the basis using QR decomposition, or are there maybe other methods in CasADi to get this basis?

Dominique

[Joel Andersson]

I am not sure I understand your question. Could you give a small example?

Be aware of the limitations of the QR factorization in CasADi - it's a very simple algorithm doing Gram-Schmidt orthogonalization without partial pivoting. 

Just discovered that it's not doing return-by-value as it should be in Python, but rather return-by-value. I've created a ticket for that:

https://github.com/casadi/casadi/issues/636

[Joel Andersson]

I meant that it is doing return-by-reference now.

[Dominique Vercammen]

Yes, I also discovered the returning by reference, also looked a little strange to me :)

Here's a small example:

A = NP.array([1,-1,-1])

As = SXMatrix(A)

Q = SXMatrix.zeros(1)

R = SXMatrix.zeros(1)

qr(As,Q,R)

print Q,R

>> [0.57735,-0.57735,-0.57735] 1.73205

So I want to find the QR decomposition of a 3-by-1 matrix A, which should return a 3-by-3 Q and a 3-by-1 R. The right part of Q, i.e., the last 2 columns, are then a basis for the null space of A transpose. The qr method however, returns a 3-by-1 Q which is only the first column.

If you would do the same in Matlab, you get this:

>> [Q,R] = qr([1 -1 -1]')

Q =

   -0.5774    0.5774    0.5774

    0.5774    0.7887   -0.2113

    0.5774   -0.2113    0.7887

R =

   -1.7321

         0

         0

with the right part of Q the null space of [1 -1 -1]:

>> null([1 -1 -1])

ans =

    0.5774    0.5774

    0.7887   -0.2113

   -0.2113    0.7887

[Joel Andersson]

I see, the current implementation assumes that the matrix is square. I have created a ticket for this:

https://github.com/casadi/casadi/issues/637

Best,

Joel

[Dominique Vercammen]

So at this point there is no way to symbolically find a null space basis using CasADi?

[Joel Andersson]

I had a look at the implementation and it looks like it does support rectangular matrices as long as there are at least as many rows as there are columns.

I don't think it's hard to extend it to calculate the rest of the columns of Q. You'll find the implementation of the function in "../../symbolic/matrix/matrix_tools.hpp".

Do you have a pointer to the algorithm?

Joel

[Dominique Vercammen]

Ok I'll take a look at the qr method.

The most used method however to find the null space basis is by doing an svd decomposition, but i don't think that is implemented in casadi?

[Joel Andersson]

No, SVD is not implemented, although I think we might do that at some point. If I remember correctly, there is a nice chain-rule for the SVD so you can actually do AD on it. Not sure if we will implement it for the scalar atomic type (SXMatrix) or only for the matrix atomic type (MX) though.

Best,

Joel

[Joel Andersson]

The qr function is now return by value again after https://github.com/casadi/casadi/issues/636 was fixed, e.g.

A = ssym("A",3,3)

Q,R = qr(A)

Joel

[Pavel Otta]

Hi Joel, I have the same problem in Maltab. Try

import casadi.*

A = SX.zeros(1,3)
A(1,1) = 1
[Q,R] = qr(A')

It results with Q (1x3), R (1x1) instead of Q (3x3), R (3x1)

Dne pátek 12. dubna 2013 14:31:24 UTC+2 Joel Andersson napsal(a):

[Pavel Otta]

Oh, I see....there are two types of QR -- 1) full and 2) thin -- and CasADi provides the second one only, right? If so, a code of full QR using Givens rotation is appended (inputs and outputs are assumed to be SX.sym).

% Golub, Van Loan: Matrix Computations, 4th edition
% QR Factorization using Givens rotation (Algorithm 5.2.4)
function [Q,R] = qrgivens(A)
import casadi.*

[m,n] = size(A);
Q = SX.eye(m);
R = SX(A);

for j = 1:n
    for i = m:-1:j+1
        [c,s] = givensrotation( R(j,j),R(i,j) );
        G = [c s; -s c];
        R([j i],j:n) = G'*R([j i],j:n);
        R(i,j) = 0; % rewrite small number by zero
        Q(j:n,[j i]) = Q(j:n,[j i])*G;
    end
end

end

% Givens rotation (Algorithm 5.1.3)
function [c,s] = givensrotation(a,b)
% Givens rotation
import casadi.*

r = if_else(abs(b) > abs(a), -a/b, -b/a);
c = if_else(b == 0, 1, if_else(abs(b) > abs(a), (1/sqrt(1+r^2))*r, 1/sqrt(1+r^2)));
s = if_else(b == 0, 0, if_else(abs(b) > abs(a), 1/sqrt(1+r^2), (1/sqrt(1+r^2))*r));

end

[Joel Andersson]

Hi,

I'm not sure if the "qr" function has been properly tested for non-square matrices.

Anyway, it's a very trivial function, the implementation is just a few lines:

https://github.com/casadi/casadi/blob/master/casadi/core/matrix.cpp#L1829

You should be able to see what's going on there by looking at the code.

Joel