How to get around sparse row deletion for least squares calculation

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?

Answer 1

[Matt J]

No, probably not. You might try pagemldivide instead, thereby avoiding sparse matrices altogether.

[Matt J]

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

Answer 2

[Bruno Luong]

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

[Matt J]

I can't see how that will be faster, though. The operation w.*A will change the sparsity of the LHS matrix, which will essentially force it to be rebuild.

[Bruno Luong]

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]

@Matt J "w.*A will change the sparsity of the LHS matrix"

No