Memory Issue - Large Poisson Matrices

Question

I'm doing some signal processing and I need to generate a Poisson matrix but the data I am working with are large enough that matlab runs out of memory the way I am currently doing it.

I've been mucking around with this for a bit and am kinda stumped at the moment, so I was hoping you guys could help me out to find a more efficient method of generating the matrix, i.e. that solves the problem of matlab running out of memory.

Anyway here is what I have so far

n = 20;

B = zeros(n^2,n^2);
for i = 1:n^2
for j = 1:n^2
    if i == j
        B(i,j) = 4;
    elseif i == j+1
        B(i,j) = -1;
    elseif i == j-1
        B(i,j) = -1;
    elseif i == j+n
        B(i,j) = -1;
    elseif i == j-n
        B(i,j) = -1;
    end
end
end

for i = n:n:n^2-1
    for j = n:n:n^2-1
        B(i+1,j) = 0;
        B(i,j+1) = 0;
    end
end

C = sparse(B);

Thanks in advance!

Answer 1

Well A+ for effort! But there is a way more efficient method for doing this. You need sparse matrices. Try something like

n = 20;
e = ones(n^2,1);
o = e;
for i = n:n:n^2-1
    o(i) = 0;
end
A = spdiags([-e -o 4*e -o -e], [-(n+1) -1 0 1 (n+1)], n^2, n^2);

if you really want to generate it yourself. I post this method so you can learn how to make banded sparse matrices using general pracitce. But for a Poisson matrix you can simply use the built in MATLAB one:

B = gallery('poisson',n);

To illustrate why you need sparse matrices, try checking the sparsity for various values of n with

sparsity = nnz(B)/prod(size(B));

Increasing the number n to somewhere around 20 is when poisson matrices really make a difference as they really are sparse (close to 1% are nonzeros). All these zeros in MATLAB are wasted space. So when you are generating B by your full for-loop method you are racking up memory! To see the difference try on your original code something like

sB = whos('B');
sC = whos('C');
disp(sA.bytes);
disp(sB.bytes);

to see that for n = 20 you get B = 1280000 bytes and C = 33928 bytes. But using this prescribed method then A = 33896 bytes! That's a difference of nearly 1.2 MB!