How to get around sparse row deletion for least squares calculation
조회 수: 1 (최근 30일)
이전 댓글 표시
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?
댓글 수: 0
답변 (2개)
Matt J
2023년 3월 6일
편집: Matt J
2023년 3월 6일
댓글 수: 1
Matt J
2023년 3월 6일
편집: Matt J
2023년 3월 6일
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(:) )
timeit( @() pagemldivide(A,b) )
Bruno Luong
2023년 3월 6일
편집: Bruno Luong
2023년 3월 6일
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
% solution after removed 10th row
xd = A(1:m-1,:)\b(1:m-1,:)
% weigthed least-square solution w
w = ones(m,1);
w (m) = 0;
xw = (w.*A)\(w.*b)
댓글 수: 3
Bruno Luong
2023년 3월 6일
편집: Bruno Luong
2023년 3월 6일
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).
참고 항목
카테고리
Help Center 및 File Exchange에서 Loops and Conditional Statements에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!