How to get around sparse row deletion for least squares calculation

1 Ansicht (letzte 30 Tage)
David Gillcrist
David Gillcrist am 6 Mär. 2023
Kommentiert: Bruno Luong am 6 Mär. 2023
I have an alrogithm that repeats a reduced basis calculation for a sparse block digaonal matrix Psi. Here Psi is composed of d number of equally sized matrices . I perform a least squares calculation for this matrix with a vector myVec to get the coefficient vector c. I then perforom an algorithm to delete a row from each of in the matrix Psi (as well as from entries in myvec). The problem I'm running into is that this is slow as I'm essentially asking for a new sparse block diagonal matrix of for each iteration, to the point where doing the least squares calculation is faster using a full matrix for Psi. Is there potentially some clever trick I can employ that allows me to delete rows from the sparse matrix without actually changing the size?

Antworten (2)

Matt J
Matt J am 6 Mär. 2023
Bearbeitet: Matt J am 6 Mär. 2023
No, probably not. You might try pagemldivide instead, thereby avoiding sparse matrices altogether.
  1 Kommentar
Matt J
Matt J am 6 Mär. 2023
Bearbeitet: Matt J am 6 Mär. 2023
With the speed-up from pagemldivide, you might be able to offset the cost of rebuilding the matrices every time:
Psi=repmat({rand(30)},1,500); blocks=cellfun(@sparse,Psi,'uni',0);
A=cat(3,Psi{:}); b=A(:,1,:); %paged form (full)
S=blkdiag(blocks{:}); %block diagonal form (sparse)
timeit( @()S\b(:) )
ans = 0.0087
timeit( @() pagemldivide(A,b) )
ans = 0.0015

Melden Sie sich an, um zu kommentieren.


Bruno Luong
Bruno Luong am 6 Mär. 2023
Bearbeitet: Bruno Luong am 6 Mär. 2023
I'm nit sure why you formalize as block sparse, since it is like solving d independent linear systems of the same size, and can be performed by pagemrdivide as Matt has pointed out.
Now back too your question, you could probably reformulate the row-deletion to a weigted lsq system:
m = 10;
n = 3;
p = 1;
format long
% Full solution
A = rand(m,n);
b = rand(m,p);
x = A\b
x = 3×1
0.092568593269531 0.674057731207846 0.304093420391084
% solution after removed 10th row
xd = A(1:m-1,:)\b(1:m-1,:)
xd = 3×1
-0.215079045247037 0.946905459896963 0.286462609971429
% weigthed least-square solution w
w = ones(m,1);
w (m) = 0;
xw = (w.*A)\(w.*b)
xw = 3×1
-0.215079045247037 0.946905459896963 0.286462609971429
  3 Kommentare
Bruno Luong
Bruno Luong am 6 Mär. 2023
Bearbeitet: Bruno Luong am 6 Mär. 2023
OP asks to keep the original SPARSE matrix, so I propose a solution for him to avoid : "I'm essentially asking for a new sparse block diagonal matrix of for each iteration".
Furthermore he could use sparse algorithms that requires matrix product vector so instead of doing (w.*A)*x, it only needs w.*(A*x).
Bruno Luong
Bruno Luong am 6 Mär. 2023
@Matt J "w.*A will change the sparsity of the LHS matrix"
No

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Loops and Conditional Statements finden Sie in Help Center und File Exchange

Produkte


Version

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by