Speeding up lsqlin to find the base of a matrix

Hello,

I have to solve the following problem over and over again for slightly different values of mu.

n = lsqlin(mu-eye(J),zeros(J,1),[],[],ones(1,J),1,[],[],[]);

The matrix mu has columns sum up to 1, values are between zero and 1 and the diagonal elements are typically above 0.9. There are only very few 0 elements in mu (although many of them are 'close' to 0, e.g. 1e-4). J is equal to 400.

You can do the same calculation using

n = null(mu-eye(J),1e-10);

n = n/sum(n);

but that's not any faster. Are there any ideas on how to speed that up? Since the solution n does not change much in the various calls, I thought about providing an initial guess, but the interior algorithm does not accept any initial guesses.

Here's a sample matrix (with J=5)

0.980472260884484 0.0169062020634941 1.31828882462499e-05 0.00712276192859210 0.00734667253541008

0.0122127547484456 0.972782440495672 1.15814051875989e-05 0.00808693210831310 0.00830831290909974

3.60311943991737e-08 4.68217920214649e-08 0.999953080612125 1.05940998862873e-07 3.71871063642989e-05

0.00445881693357994 0.00660554154798086 1.18490160022582e-05 0.976444331730883 0.0148592752450766

0.00285613140229620 0.00370576907106219 1.03060784389177e-05 0.00834586829121426 0.969448552204049

Thanks!

[Matt J]

I have to solve the following problem over and over again for slightly different values of mu.

Are the different mu available at the same time?

[John D'Errico]

There are at least a couple of other ways I could describe to solve for the null space of a matrix. It would be easier if you would provide a sample matrix to play with, and compare methods, without needing to describe in detail how to solve it for each method.

[Tintin Milou]

The different mu are not available at the same time. I'll provide a sample matrix to play with.

Answer 1

[Bruno Luong]

Can you try this:

A = mu-eye(size(A));
[Q,R,p] = qr(A,'vector');
n = [R(1:end-1,1:end-1)\R(1:end-1,end); -1];
n(p) = n/sum(n)

[Tintin Milou]

That's very good. Thanks! It cuts the time by about 2/3.

[Bruno Luong]

It seems faster

A=rand(400,399)*rand(399,400);
ntests = 10;

tic;
for k=1:ntests
    n=null(A); 
    n = n/sum(n);
end
toc

Elapsed time is 0.261532 seconds.

tic;
for k=1:ntests
    [Q,R,p] = qr(A,'vector');
    x = [R(1:end-1,1:end-1)\R(1:end-1,end); -1];
    x(p) = x/sum(x);
end
toc

Elapsed time is 0.079755 seconds.

norm(n-x)/norm(n)

ans = 3.4488e-08

[John D'Errico]

QR would have been my best guess. It is less computationally less intensive than an SVD, so it will be faster there. RREF is another possibilty to yield a nullspace, but RREF is written in MATLAB as I recall, not compiled code. So RREF is relatvely slow.

[Tintin Milou]

Update:

I'm using this piece of code within a call to lsqnonlin. That's why it is evaluated so many times.

The three methods (lsqlin, null, qr) give slightly different values for n, as also reported by Bruno. Due to the differences, the three methods yield a different solution in the outer maximization problem. In particular, with qr, the outer solver does not always converge, while it does for the other two. When it converges, it is faster, but I'm not sure how I can make sure that it always converges.

It's puzzling to me that such small differences between the two methods yield very different results for the outer maximization problem. I know that I can manage the precision when using lsqlin (TolFun is set to 1e-12). Is there a similar thing for qr?

[Bruno Luong]

did you get any MATLAB warning "Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = xxxxxx." at this line

n = [R(1:end-1,1:end-1)\R(1:end-1,end); -1]

[Bruno Luong]

This is probably less prone to issue with dim(kernel(A))>1. I hope it is still fast.

A = mu-eye(size(A));

[Q,R,p]=qr(A,0);
n=size(R,1)+1;
d=abs(R(1:n:end-n));
tol=1e-10;
j=find(d>=tol*d(1),1,'last');
x=[R(1:j,1:j)\R(1:j,end);zeros(n-j-2,1);-1];
x(p) = x/norm(x)

[Bruno Luong]

@Tintin Milou "I'm using this piece of code within a call to lsqnonlin. That's why it is evaluated so many times."

How do you deal the fact that n can flip the sign and it could in turn create a discontinuity that lsqnonlin might not like?

[Tintin Milou]

No, I did not get any warnings.

n itself is a vector with all values between 0 and 1, and the sum of their elements equal to 1. So it never flips signs.

Your new method gives me the same results as the old one.

I'm re-adjusting my outer loop and keep you posted. Thanks for your help.

[Tintin Milou]

Update: Before calling fsolve for the outer loop, I now add an fsolve on a small-scale problem to get a reasonable initial guess for the complete problem. With this improved initial guess, the code runs more smoothly and the initial suggestion by Bruno Luong is the best solution I have found. Thanks again!

[Bruno Luong]

Thanks for the update. It is puzzled me that the outer loop is that sensitive to numerical error.