Is it possible to save only half of a symmetric matrix to save the memory?
There is a large matrix that is used in Ax=b
type problem. A is symmetric. Is there any algorithm to let us save only half of the matrix a开发者_开发知识库nd do operation like x=A\b
on it?
You'll only save half the memory, but you can do this by creating a flat version of the matrix, saving that, then unflattening it. The extra time required probably doesn't make the saving worthwhile, mind:
% pretend this is symettric...
A = rand(10, 10);
% store it as a flat list
flatA = [];
for k = 1:size(A,1)
flatA = [flatA; A([1:k],k)];
end
save('flatA.mat', 'flatA');
% undo
N = (-1 + (1+4*2*length(flatA))^0.5)/2;
newA = NaN(N,N);
cur = 1;
for k = 1:N
len = k;
newA(1:len, k) = flatA(cur:cur+len-1);
cur = cur + len;
end
% real A cannot have NaNs or this trick fails
newA(isnan(newA)) = newA(isnan(newA'))';
Here's an idea, but I haven't tested it out. If your symmetric matrix is positive definite, do a Cholesky decomposition of the symmetric matrix, A, giving you A = U*U'. If you store U as a sparse matrix using MATLAB's builtin sparse matrix, you have everything you need, and you've used roughly half the memory. Since its using MATLAB's sparse matrix type, you have operate on it using standard MATLAB functions, as long as you remember that A = U*U'
For example, to compute A*x = b, use x = U'\U\b. Unlike the other proposed solutions, MATLAB will never be actually using a full matrix internally, and will even use accelerated solvers that will take advantage of the fact that you're only solving with triangular. The cost is that to solve a single system, you've actually running the backslash operator twice (see above). However, that's the price you pay for never instantiating the full matrix.
If you extract the upper triangular part and convert to a sparse matrix, it should save memory. This technique is reasonably quick.
% Create a symmetric matrix
A = rand(1000);
A = A + A';
% Extract the upper triangular part
B = sparse(triu(A)) % This is the important line, the rest is just checking.
% Undo
BB = full(B);
C = BB + BB' - diag(diag(BB));
% Check correct
assert(all(A(:) == C(:)));
% Save matrices
save('A.mat', 'A');
save('B.mat', 'B');
% Get file sizes
infoA = dir('A.mat'); infoA.bytes
infoB = dir('B.mat'); infoB.bytes
EDIT to clarify things for woodchips
The above solution demonstrates a way of saving matrices with less file space. The matrix B
also takes up less memory than the original matrix A
. If you want to do linear algebraic operations on B
, that works perfectly well. Compare
b = rand(1000);
fullTriUA = triu(A);
sparseTriUA = sparse(fullTriUA); %same as B above
fullResult = fullTriUA\b;
sparseResult = sparseTriUA\b;
assert(all(fullResult(:) == sparseResult(:)))
精彩评论