Out of memory using mrdivide

8 views (last 30 days)
Bas
Bas on 1 Jun 2011
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
  1 Comment
Jan
Jan on 1 Jun 2011
Do you have 1 element per line? Then A/B = (B'\A')' could be done manually also.

Sign in to comment.

Accepted Answer

Ivan van der Kroon
Ivan van der Kroon on 1 Jun 2011
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.
  1 Comment
Bas
Bas on 4 Jun 2011
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

Sign in to comment.

More Answers (6)

Ivan van der Kroon
Ivan van der Kroon on 5 Jun 2011
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.
  3 Comments
Ivan van der Kroon
Ivan van der Kroon on 6 Jun 2011
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
Bas on 6 Jun 2011
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.

Sign in to comment.


Sean de Wolski
Sean de Wolski on 1 Jun 2011
Buy more RAM and upgrade to a 64bit OS.
  1 Comment
Bas
Bas on 6 Jun 2011
I have a 64bit OS and 4GB of RAM.

Sign in to comment.


Image Analyst
Image Analyst on 4 Jun 2011
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;
  1 Comment
Bas
Bas on 4 Jun 2011
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.

Sign in to comment.


Ivan van der Kroon
Ivan van der Kroon on 7 Jun 2011
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.
  10 Comments
Bas
Bas on 13 Jun 2011
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
Ivan van der Kroon on 15 Jun 2011
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?

Sign in to comment.


Bjorn Gustavsson
Bjorn Gustavsson on 7 Jun 2011
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,
  3 Comments
Bjorn Gustavsson
Bjorn Gustavsson on 9 Jun 2011
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
Bas on 13 Jun 2011
Bjorn, Do you maybe have a reference about your solution to an inhomogeneous heat equation? It seems promissing.

Sign in to comment.


Bjorn Gustavsson
Bjorn Gustavsson on 13 Jun 2011
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!