Speeding up lsqlin to find the base of a matrix

2 visualizaciones (últimos 30 días)
Tintin Milou
Tintin Milou el 4 de Ag. de 2022
Comentada: Bruno Luong el 7 de Ag. de 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 comentarios
John D'Errico
John D'Errico el 4 de Ag. de 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 el 4 de Ag. de 2022
The different mu are not available at the same time. I'll provide a sample matrix to play with.

Iniciar sesión para comentar.

Respuesta aceptada

Bruno Luong
Bruno Luong el 4 de Ag. de 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 comentarios
Tintin Milou
Tintin Milou el 7 de Ag. de 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 el 7 de Ag. de 2022
Thanks for the update. It is puzzled me that the outer loop is that sensitive to numerical error.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Parallel for-Loops (parfor) en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by