"Out of Memory" when inverting a symbolic matrix

I need to set up a program using symbolic matrices. However, the mathematics apparantly becomes complicated enough that taking the inverse of a 3x3 symbolic matrix takes up too much memory. I wanted to ask how I can circumvent the memory issue. Following is the relevant code segment leading to the problem:

v=['[' sprintf('conca%d,',1:51)];
v(end)=']';
Cconc1 = str2sym(v);
clear v
v=['[' sprintf('concb%d,',1:51)];
v(end)=']';
Cconc2 = str2sym(v);
clear v
v=['[' sprintf('kk%d,',1:2)];
v(end)=']';
KK = str2sym(v);
KK = KK';
Cconc1 = Cconc1';
Cconc2 = Cconc2';

% Setting up symbolic CC matrices
    CC  = [Cconc1,  KK(1)*Cconc1.*Cconc2,  KK(1)*KK(2)*Cconc1.*Cconc2.^2];
    
% Concerning EE calculation:
%   No difference between using \ or inv() on memory usage
%   Splitting the equation into smaller sections as below also does not assist

preEEinv1 = (CC'*CC);
preEEinv2 = (inv(preEEinv1));
clear preEEinv1
EE = (preEEinv2*(CC'*DD));
% Note - CC is a 51x3 matrix, KK is a 2x1 vector, DD is a 51x1 vector

% previous attempts -->     1)   EE = (CC'*CC)\(CC'*DD);
%                           2)   EE = inv(CC'*CC)*(CC'*DD);

The following is the error message returned when executing the code

Error using symengine
Out of memory.

Error in sym/privUnaryOp (line 1045)
            Csym = mupadmex(op,args{1}.s,varargin{:});

Error in sym/inv (line 22)
X = privUnaryOp(A, 'symobj::inv');

Error in sym_r (line 98)
preEEinv2 = (inv(preEEinv1));

As I am not used to getting out of memory errors on my PC, I'm a little out of my depth. Also, if you require information regarding my hardware, I am working with 16 GB RAM, a i7 3.60 GHz processor, and 65 GB virtual space.

The error line indication is always where the inverse is taken. I don't make unnecessary copies of data, and cannot use sparse matrices. The tall function does not work for the symbolic toolbox either. I have tried inverting much simpler symbolic matrices to make sure the function is sound, and it worked. Is there a way I can bypass or solve the memory issue?

Answer 1

[Walter Roberson]

syms EE [3 3]
EEinv = inv(EE);
preEEinv2 = subs(EEinv, EE, preEEinv1);

That is, instead of directly taking the inverse of your calculated 3 x 3 matrix, instead take the inverse of a sample 3 x 3 matrix, and once you have that, substitute in the actual matrix entries for the place-holders. This is typically much faster.

By the way: it appears to me that you are probably using real-valued quantities. If so, then I recommend that you add assumptions, such as

assume(KK,'real')
assume(Cconc1 >= 0)

as appropriate for your situation. Your preEEinv1 has a lot of conj() calls that can be removed if you are working with real values.

[Wikus Bergh]

Thank you for your answer Walter, I appreciate the quick response. I tested it on dummy matrices, and it works perfectly. You made my life much easier now, by giving me a method I can use in multiple code segments to increase performance.

[Walter Roberson]

One further benefit is that if you have a series of inverses of 2 x 2, 3 x 3, or 4 x 4 to do, then you can pre-compute the inv() of the template, and then just have to subs() afterwards.

I do not, however, recommend doing inv() of a symbolic 5 x 5 or later. 5 x 5 would be the absolute upper limit, for the case where you were willing to give it several hours to pre-compute the inverse (but then the next step with the subs() will probably give you something that is too large to really work with.)

[Wikus Bergh]

Hi Walter, out of curiosity, what should be done when you need to invert a larger symbolic matrix? Is is only then dependent on how much time you have, or can the problem be broken down into smaller pieces and analysed as such?

[Walter Roberson]

Correction: up to 6 x 6 is fast. I am not certain if 7 x 7 is fast or not (I made a mistake and need to redo it). 8 x 8 is slow.

[Walter Roberson]

https://en.wikipedia.org/wiki/Invertible_matrix#In_relation_to_its_adjugate

inv(S) = 1/det(S) * adj(S)

where

https://en.wikipedia.org/wiki/Adjugate_matrix#Definition

adj(S) = C(S)^T where C(S) is the cofactor matrix and ^T is transpose.

So you can build up the inverse in pieces by first calculating the determinanent and then calculating the elements of the cofactor matrix; you can do the division for each element as you go, perform appropriate substitutions, and save the result in a separate file if you need to.

But if you do not have enough memory to store the entire inverse, then what are you going to do with the inverse? You would not typically calculate the inverse of a matrix just for fun: you would normally be wanting to do a calculation with it, multiplying it by another matrix -- and to do that, you would need a row (or column) of the inverse for each different output of the multiplied matrix.

Basically you are always on the edge of disaster for that kind of calculation: if an N x N matrix strains your resources, then an (N+1) x (N+1) matrix or (N+2) x (N+2) matrix would overwhelm it.

I guess there might be some approaches involving multiple workers in parallel, and shared memory (or shared file) system, and a system of locks to ensure that calculations are not re-done... it would probably have to be SNMP. Maybe not a system of locks, but instead a system of queues... that might make more sense.