Speeding up lsqlin to find the base of a matrix

2 Ansichten (letzte 30 Tage)
Tintin Milou
Tintin Milou am 4 Aug. 2022
Kommentiert: Bruno Luong am 7 Aug. 2022
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!
  3 Kommentare
John D'Errico
John D'Errico am 4 Aug. 2022
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
Tintin Milou am 4 Aug. 2022
The different mu are not available at the same time. I'll provide a sample matrix to play with.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Bruno Luong
Bruno Luong am 4 Aug. 2022
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)
  10 Kommentare
Tintin Milou
Tintin Milou am 7 Aug. 2022
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
Bruno Luong am 7 Aug. 2022
Thanks for the update. It is puzzled me that the outer loop is that sensitive to numerical error.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Parallel for-Loops (parfor) finden Sie in Help Center und File Exchange

Tags

Produkte


Version

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by