Out of memory using mrdivide

2 次查看（过去 30 天）
Bas
Bas 2011-6-1
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 个评论
Jan
Jan 2011-6-1
Do you have 1 element per line? Then A/B = (B'\A')' could be done manually also.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

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.
  1 个评论
Bas
Bas 2011-6-4
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

请先登录，再进行评论。

更多回答（6 个）

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.
  3 个评论
显示 1更早的评论
Ivan van der Kroon
Ivan van der Kroon 2011-6-6
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 2011-6-6
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
Sean de Wolski 2011-6-1
Buy more RAM and upgrade to a 64bit OS.
Image Analyst
Image Analyst 2011-6-4
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 个评论
Bas
Bas 2011-6-4
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
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 个评论
显示 8更早的评论
Bas
Bas 2011-6-13
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 2011-6-15
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
Bjorn Gustavsson 2011-6-7
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 个评论
显示 1更早的评论
Bjorn Gustavsson
Bjorn Gustavsson 2011-6-9
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 2011-6-13
Bjorn, Do you maybe have a reference about your solution to an inhomogeneous heat equation? It seems promissing.

请先登录，再进行评论。

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

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Analysis 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by