Thread Subject: factorize singular symmetric matrix

Let A a real symmetric positive - but *not* definite - matrix.

What is the good way to factorize P'*A*P where P is (any - but to be found) orthogonal basis of Im(A)? In other word, am I obligate to perform QR factorization of A and throw out the "smaller" vectors. This method cannot take any advantage of the symmetric positiveness of A.

My A might be sparse. So I would also prefer to avoid factorization filling, but let alone this constraint for the moment.

Bruno

maybe something like this:

% this tolerance is used to determine the rank of the matrix
tol=1e-12;

% find the rank by means of Cholesky factorization
[L,p]=chol(A);
p=p-1;
i=0;
while sum(abs(L(p-i,end-i:end)).^2)<tol
   i=i+1;
end
rk=p-i;

% find the orthogonal vectors basis for the space generated by A
U=zeros(M,rk);
for i=1:rk
   temp=A(:,i);
   for j=1:i-1
      temp=temp-U(:,j)*(U(:,j)'*temp);
   end
   U(:,i)=temp/sqrt(sum(abs(temp).^2));
end


Stefano


"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h9sqsb$7gk$1@fred.mathworks.com>...
> Let A a real symmetric positive - but *not* definite - matrix.
>
> What is the good way to factorize P'*A*P where P is (any - but to be found) orthogonal basis of Im(A)? In other word, am I obligate to perform QR factorization of A and throw out the "smaller" vectors. This method cannot take any advantage of the symmetric positiveness of A.
>
> My A might be sparse. So I would also prefer to avoid factorization filling, but let alone this constraint for the moment.
>
> Bruno

"Stefano " <s.mangione@gmail.com> wrote in message <h9t8h8$3tp$1@fred.mathworks.com>...

>
> % find the rank by means of Cholesky factorization
> [L,p]=chol(A);

Do you think cholesky can overcome matrix such like

A=[0 0;
     0 1]

Bruno

I assumed that there were no zeros on the diagonal... substitute that line with:

[L,p]=chol(A+tol^2*eye(M));

Stefano

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h9t9bf$o6b$1@fred.mathworks.com>...
> "Stefano " <s.mangione@gmail.com> wrote in message <h9t8h8$3tp$1@fred.mathworks.com>...
>
> >
> > % find the rank by means of Cholesky factorization
> > [L,p]=chol(A);
>
> Do you think cholesky can overcome matrix such like
>
> A=[0 0;
> 0 1]
>
> Bruno

"Stefano " <s.mangione@gmail.com> wrote in message <h9tabe$2e$1@fred.mathworks.com>...
> I assumed that there were no zeros on the diagonal... substitute that line with:
>
> [L,p]=chol(A+tol^2*eye(M));
>

The problem is not the diagonal, the problem is cholesky is unable to estimate in reliable way the rank of the matrix. Another example:

 A=[1 4 4;
      4 20 20;
      4 20 20]

This matrix has rank = 2.

Bruno

Sorry bad example previously, please use this one:

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

[L,p]=chol(A)

rank(A) % <- 2

Bruno

shame on me, I don't know how to get a good estimate of the rank :(

moreover, the Gram Schmidt implementation above is wrong since it adds a vector to the basis even if it is zero-norm.

Stefano

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h9tbqo$q6v$1@fred.mathworks.com>...
> Sorry bad example previously, please use this one:
>
> A = [1 2 0 0;
> 2 4 0 0;
> 0 0 1 2;
> 0 0 2 4]
>
> [L,p]=chol(A)
>
> rank(A) % <- 2
>
> Bruno

"Stefano " <s.mangione@gmail.com> wrote in message <h9tgq4$2aa$1@fred.mathworks.com>...
> shame on me, I don't know how to get a good estimate of the rank :(
>

Don't be ashamed, the problem is harder than it seems to be (at least for me).

Bruno

at last, I think I've found a way to get with chol() the same reliability of the rank() function (seems faster, tested with the script at http://tinyurl.com/yevnhmk ):

tol=eps*norm(A)*numel(A);
[L,ignored]=chol(A+tol*eye(M));
rk=sum(abs(diag(L))>sqrt(sqrt(tol)));

in order to get the orthogonal basis, though, I wanted to use the Gram-Schmidt algorithm; now, when the matrix is not block-diagonal, the algorithm results indeed faster than a QR; it turns out that when the matrix is block-diagonal the loop is slower than using the more straightforward QR factorization. In order to get a faster algorithm, a fast way to test for block-diagonality is needed; then it should be trivial to get a faster algorithm.

Stefano

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h9thjr$onk$1@fred.mathworks.com>...
> "Stefano " <s.mangione@gmail.com> wrote in message <h9tgq4$2aa$1@fred.mathworks.com>...
> > shame on me, I don't know how to get a good estimate of the rank :(
> >
>
> Don't be ashamed, the problem is harder than it seems to be (at least for me).
>
> Bruno

Stefano,

Note that Gram-Schmidth is nothing more than an unstable way of doing QR factorization (the later uses Householder reflection or Given rotation to "orthogonalize" the vectors).

Bruno

Stefano, I find this paper where one can find deep discussion about rank estimation based on Cholesky and QR. http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf

Time for some reading.

Bruno

Another interesting paper.
http://www.emis.de/journals/ETNA/vol.17.2004/pp76-92.dir/pp76-92.pdf

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h9tpbc$7qh$1@fred.mathworks.com>...
> Stefano, I find this paper where one can find deep discussion about rank estimation based on Cholesky and QR. http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf

A quick read through the paper, the author seems to suggest that Cholesky with or without pivoting are not a reliable method to estimate the rank. I try one of his example (Example 2.1) without pivoting (BTW do we have cholesky with pivoting in Matlab?) and it might overestimate the rank. He recommended QR as a better method.

Bruno

Thank you for the links. Still I think that my implementation is as robust as matlab's rank(), while being faster (and applicable only to symmetric hermitian matrices)

Stefano

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h9u482$a9l$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h9tpbc$7qh$1@fred.mathworks.com>...
> > Stefano, I find this paper where one can find deep discussion about rank estimation based on Cholesky and QR. http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf
>
> A quick read through the paper, the author seems to suggest that Cholesky with or without pivoting are not a reliable method to estimate the rank. I try one of his example (Example 2.1) without pivoting (BTW do we have cholesky with pivoting in Matlab?) and it might overestimate the rank. He recommended QR as a better method.
>
> Bruno