Matlab: construction of very large sparse band matrix with as little memory and computation

Question

I need to construct a large NxN sparse band matrix A with N = 570*720 = 410400 (# of image pixels).

Mathematically, A(m,n) = C1 * exp(-|m-n|^2); m = 1:N, n = 1:N

Basically its a Gaussian function evaluated at each row with row index being the mean and some arbitrary but small standard deviation.

For N = 100, it looks like:

Unfortunately, it runs very slow for N = 410400 due to unnecessary computations.

1) using for loop

A = sparse(N,N);
for i=1:N
A(i,:) = normpdf(1:N, i, 30);
end

This is wasteful and slow due to calling normpdf N times.

2) compute normpdf once for 1:2N with mean at N and then circularly shift the row in A with appropriate mean based on index. circshift in matlab can't shift a matrix column wise with different shift sizes. Again will need to call circshift N times --> wasteful.

3) use mvnpdf perhaps but it needs input vectors and generating these with meshgrids will
consume lot of memory.

4) use bsxfun with user defined gaussian function (with fixed SD) accepting two parameters as bsxfun does not take >3 arguments.

Any ideas on how this can achieved efficiently?

Thanks

Answer 1

May I suggest another approach altogether.

What would be against using a 2*N-by-1 vector, that can be indexed with a simple transformation on the indices? Like so:

% Oli's solution, and your request: NxN matrix
N   = 100;
s   = 30;    
pdf = @(x) 1/(sqrt(2*pi)*s)*exp(-0.5*(bsxfun(@minus,x(:),1:N)/s).^2);
A   = pdf(1:N);

% My solution: 2*N x 1 vector
B = exp(-0.5*((-N:N)/s).^2) / s/sqrt(2*pi);

The trick is to find a nice general indexing rule. Here's how to do it:

% Indexing goes like this:
fromB = @(ii,jj) B(N+1 + bsxfun(@minus, jj(:), ii)).';

ii = 30; 
jj = 23;

from_A = A(ii,jj)
from_B = fromB(ii,jj)

ii = 1:2;    
jj = 4:6;      

from_A = A(ii, jj)
from_B = fromB(ii,jj)

Results:

from_A =
   0.012940955690785
from_B =
   0.012940955690785

from_A =
   0.013231751582567   0.013180394696194   0.013114657203398
   0.013268557543798   0.013231751582567   0.013180394696194
from_B =
   0.013231751582567   0.013180394696194   0.013114657203398
   0.013268557543798   0.013231751582567   0.013180394696194

I'm sure you can figure out how to do things like colon-indexing and using the end keyword :)

Answer 2

First of all, you really do not need a sparse matrix unless you have at least 50% of zeros but your matrix is full.

Consider the pdf of a normal

You can easily implement it including a call to bsxfun():

N   = 100;
s   = 30;
m   = 1:N;
pdf = @(x) 1/(sqrt(2*pi)*s)*exp(-0.5*(bsxfun(@minus,x(:),m)/s).^2);
B   = pdf(1:N);

A simple example with mean = 1 and sigma = 30 will clarify:

pdf = @(x) 1/(sqrt(2*pi)*30)*exp(-0.5*((x-1)/30).^2);
pdf(1)
ans =
        0.0132980760133811

normpdf(1, 1, 30)
ans =
        0.0132980760133811