Out of memory using mrdivide

Hello,

I am dividing two sparse matrices.

A / B

I get the message: ??? Error using ==> mrdivide Out of memory. Type HELP MEMORY for your options.

The Largest available Contiguous Free Blocks on my computer has 2046 MB. Both sparse matrices have 36288 non-zero elements. The size of the full matrices is 36288*36288. I do not want an element by element division. Does anyone know how to solve this problem?

Thanks!

Bas

[Jan]

Do you have 1 element per line? Then A/B = (B'\A')' could be done manually also.

Answer 1

[Ivan van der Kroon]

I agree with Jan; if you have the same number of non-zero elements you probably filled the main diagonal. In that case spdiags works great, e.g.

n=36288;
d=diag(A)./diag(B); %diagonal of the requested matrix
C=spdiags(d,0,n,n);

If they are not on the main diagonal, just be careful about which diagonal to fill.

[Bas]

Actually my problem is more complicated than I outlined. The A matrix has indeed only the main diagonal filled, the B matrix has the main diagonal filled plus one sub-(i<0) and one super(j>0) diagonal, where: i = -j.

Do you or Jan know of a trick how to divide B from A?

Thanks in advance!

Bas

Answer 2

[Ivan van der Kroon]

Some follow-ups:

a) Where do you need this for? Your result will lose much of its sparcity. Maybe if you post your entire problem, some more practical solutions can be suggested. You probably want to solve some ode?

b) About the sparcity; you can predict which entries are going to be filled. Since A is diagonal the entries are equal to the entries in the inverse of B. B has three diagonals; one at 0 (i.e. the main diagonal, as in spdiags) and two others at which position? This is very important because the inverse is now filled at the multiples of the location of these diagonals. For example, when they are at -1 and 1 the inverse is a full matrix, when it is -2 and 2 the inverse is empty for its odd diagonals (i.e. 0, -2, 2, -4, 4,… are filled), when they are at -3 and 3 the inverse diagonals are at 0, -3, 3, -6, 6,... and so on. Especially when the diagonals are close to the main diagonal, I would go into a) and stay away from this matrix AB^-1.

[Bas]

Good idea. I'll explain the full problem, maybe there are other solutions than I am thinking off.

I am simulating dispersal of insects via a Fokker Planck diffusion process. During the dispersal, insects are trapped in a grid of pitfall traps.

My Partial Differential Equation has the form:

(d/dt)N = D((d2/dx2)N + (d2/dy2)N) - (a + u)N

in which: N are the number of insects at time t and location x,y; D (m2/day) is the diffusion constant, or motility; a (1/day) is the removal rate of insects in traps which is location specific; u (1/day) is the removal rate of insects due to death and is not location specific.

D, a and u are unknown parameters which I need to estimate from the a data set of a mark-recapture experiment. We tried, but could not solve this PDE analytically.

I used the forward central (FC) method to numerically solve the PDE, but this method requires a small time step of integration and has a rather large numerical error.

The Alternating Difference Equation (ADI) method on the contrary allows a larger time step and has a smaller error. The ADE, however, can only be solved in a matrix, vector form. And there resides my problem!

I will look into your suggestion (b) for a solution. Mean while I'll keep my ears open for other ideas.

[Ivan van der Kroon]

Can you specify where this matrix division comes in? I assume size(N)=size(a)=[n,n]; size(D)=size(u)=[1,1]; n=36288;

Then it reads dN/dt=nabla^2 N - (a(x,y)+u)*N where nabla is also known as the del operator. The diffusion equation is usually solved by separation of variables. The problem here is the source term a. How 'nice' is this function?

[Bas]

You're quite right: size(N)=[n,1]; size(a)=[n,n], size(D)=size(u)=[1,1]; n=36288.

a(x,y)=(g*D)/(dx*dy), where dx and dy are the spatial discretizations; dx=dy=1. g is a scaling parameter that is used to approximate a(x,y).

We separated variables and concluded that since g is unknown we cannot analytically solve the PDE.

The matrix division comes in after transforming the PDE to matrix/vector form. This equation has two matrix divisions of which I show the first one:

Nt+dt=(I-c*I*u-c*A)/(I-c*Qy)...

here c represents different constants with size=[1,1]. I is the identity materix. The diagonal of A is filled with the value g for those elements that correspond to trap locations in the field. Qy is the transition matrix in the y-direction and has -2 on the main diagonal and 1 on the k=189 and k=-189 diagonal. In my question above

A=I-c*I*u-c*A and B=I-c*Qy.

Hope this helps.

Answer 3

[Sean de Wolski]

Buy more RAM and upgrade to a 64bit OS.

[Bas]

I have a 64bit OS and 4GB of RAM.

Answer 4

[Image Analyst]

Are you sure you want to do a matrix division (like multiplying by the matrix inverse) rather than an "element-by-element" division of A by B. If you want every element of A divided by the corresponding element of B you'd use dot slash instead of slash

elementByElementDivision = A ./ B;

[Bas]

No, it's indeed an awful division of two matrixes. For those interested: I am trying to implement an Alternating Difference Equation for a numerical solution of a diffusion process.

Answer 5

[Ivan van der Kroon]

Short about inv(B); it has 191*2+1 diagonals (36288/189=192, but this is for one of the triangles and includes the main diagonal). It is also symmetric and off course you can divide the matrix in C=192x192 sections of I=189x189 identity matrices, such that inv(B)=kron(C,eye(189)); Hence there will only be 192*192-192=36672 distinct values. The entries will be

a=-192:-1;
C=zeros(192);
C(n,n:end)=n*a(n:end);%e.g. with a for-loop and make it symmetric afterward)
C=C/193;

You can verify this for instance for 1=B(1,1)*inv(B)(1,1)+B(190,1)*inv(B)(1,190)=-2*-192/193+1*-191/193=193/193=1.

But I just don't think this is practical because you will still be out of memory when you try to perform kron, i.e. making the big matrix inv(B); However, you know all the entries now.

I don't get the a(x,y) since N is 1-D and a doesn't show up in your discrete time integration.

[Ivan van der Kroon]

Sorry, I just made a mistake; this is of course for inv(Qy), but I'm sure you can find the entries of inv(B) yourself.

[Bas]

Hello Ivan, I tried your approach for a smaller matrix, but don't find that inv(Qy) = kron(C,eye(n)). Can you verify if I got you right?

Regarding your question: a(x,y) is approximated by (g*D)/(dx*dy) for those locations where there are traps and is zero otherwise. In the discrete time integration the recapture of insects is represented by c*A*N, where c is a constant, A is a matrix with the value g on the main diagonal for those locations that contain a trap, and N is the vector that contains the number of beetles on each location.

% Here I try to verify that inv(Qy) = kron(C,eye(n))

x=18; % number of nodes in the x-direction (this was 189)

y=19; % number of nodes in the y-direction (this was 192)

Qy = sparse(x*y,x*y);

for i = 1:(x*y-x)

Qy(i,i+x)=1;

Qy(i,i)=-2;

Qy(i+x,i)=1;

end

a=-y:-1;

C=zeros(y);

for i=1:y

C(i,i:end) = i*a(i:end);

C(i:end,i) = i*a(i:end);

end

C=C/y+1;

sum(sum(inv(Qy) - kron(C,eye(x)))); % this does not yield a zero

[Bas]

I just realize that inv(Qy) does not yield memory problems. It goes wrong at inv(Qx). Qx has -2 on the k=0 diagonal. And has 1 on the k=-1 and k=1 diagonal except for the i*(x,x+1) and i*(x+1,x) elements where i=1:y-1.

I assume that inv(Qx) = kron(C,eye(n)) should also be valid for Qx.

[Ivan van der Kroon]

The matrix Qy you produce does not have a full main diagonal. You can check it with spy(Qy). It is easier to make it like this Qy = spdiags(ones(x*y,1)*[1 -2 1],[-x 0 x],x*y,x*y); It just checked it and get something like 1e-14 for norm(inv(Qy) - kron(C,eye(x)));

Can you split a(x,y)=a1(x)*a2(y)? Can you also specify what (d2/dx2)N and (d2/dy2)N are since N is just 1-D? If I interpret your equation as it is, you can just solve it with ode45. Can you you maybe specify again what exactly the unknowns are?

[Bjorn Gustavsson]

I guess N is the 2-D density of insects. So Bas wants to solve for the timevariation of N(x,y) that is get N(x,y,t)

[Ivan van der Kroon]

Well, Bas commented "size(N)=[n,1]; size(a)=[n,n], size(D)=size(u)=[1,1]; n=36288."

[Bjorn Gustavsson]

That seems odd to me, but I might be confusing things. The PDE sure looks like a 2-D diffusion equation with sources and losses...

[Bas]

Ivan, I don't see how you get 1e-14 for norm(inv(Qy) - kron(C,eye(x)));

I use the following code and get 18.1045 for norm(inv(Qy) - kron(C,eye(x)));

x=20; % number of nodes in the x-direction (this was 189)

y=20; % number of nodes in the y-direction (this was 192)

Qy = spdiags(ones(x*y,1)*[1 -2 1],[-x 0 x],x*y,x*y);

a=-y:-1;

C=zeros(y);

for i=1:y

C(i,i:end) = i*a(i:end);

C(i:end,i) = i*a(i:end);

end

C=C/y+1;

norm(inv(Qy) - kron(C,eye(x)))

[Bas]

OK, I see that the 1-D of N causes trouble. I try to explain it. Supose that my space of insects consists of 3x4 nods: size(space) = [3,4]. Then space = reshape(N,4,3). So the elements in N are [space(1,1); space(1,2); space(1,3); space(1,4); space(2,1);...space(3,4)].

To make it all clear I'll present my PDE again with the unknowns:

dN/dt=D*nabla^2*N - (a(x,y)+u)*N

a(x,y) is congruent to (triple bar): (g*D)/(dx*dy). So a(x,y) is not a nice function. a=(g*D)/(dx*dy) for the nods with traps, for the other nods a=0;

The unkowns in this PDE are D, g and u. These values are being estimated from data using minimum likelihoods.

a(x,y) can be split into a1(x) and a2(y). Then a1(x) = a2(y). Why is this important Ivan?

[Ivan van der Kroon]

You don't get the correct matrix because you do not normalize it correctly: C=C/(y+1); You divided by y and added one. For both x and y at 20 I get 1.0073e-013 on my win7x64 machine.

I assume you are able to take the correct Laplacian although you have reshaped this vector. What are the sizes of D,G,u,dx,dy? Are they available or unknowns (I don't quite understand what you mean by estimating)?

Whether a can be split makes the number of variables much smaller (for n it now is 2n instead of n^2). Next I figured that you can perform separation of variables on N, but you clarified now what N is? So, how do I interpret a as N is actually 2-D, i.e what does a looks like after reshaping?

Answer 6

[Bjorn Gustavsson]

Bas, why cant you solve the PDE analytically - or at least "analytically enough"? A diffusion is just a convolution with a Gaussian Kernel, so the solution of your diffusion equation should be something not too dissimilar. See for example: http://en.wikipedia.org/wiki/Heat_equation#Fundamental_solutions. That should be possible to use to build a discretized model of to use for parameter estimation.

HTH,

[Bas]

Hello Bjorn, Do you think a convolution also works when objects are taken out of the diffusion process at certain locations. The removal of insects due to trapping is a significant process in my case that I cannot ignore.

[Bjorn Gustavsson]

Yes, I think so. I've used a slightly modified 3-D version of the Greens function solution to "Inhomogeneous heat equation,

Problem on (-∞,∞) homogeneous initial conditions" on that wikipedia page for a time-varying parameter estimation problem. It should be possible to combine that solutions with the one for the initial value problem. I don't know how "discrete" your insect density is - my problem was for a truly continous gas. But for the PDE you wrote above it should work.

[Bas]

Bjorn, Do you maybe have a reference about your solution to an inhomogeneous heat equation? It seems promissing.

Answer 7

[Bjorn Gustavsson]

Not realy anything better than the references on the wikipedia pages. But http://eqworld.ipmnet.ru/en/solutions/lpde/heat-toc.pdf contains a few more examples. That the solutions solve the heat equation should be possible to prove by differentiating the solutions.