# Out of memory using mrdivide

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

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.
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?
Bas

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.
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.

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

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;
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.

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.
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?

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

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.