4

The problem is as follows:
I have a set of 20 equations:

r1 = s1*h1_1 + s2*h1_2 + ... s20*h1_20

r2 = s1*h2_1 + ...

...

r20 = s1*h20_1 + ...

where r, s and h are matrices and * symbolizes pointwise product. It can be rewritten in the matrix form R = H*S. I want to solve this equation for S - so I need to calculate inv(H)*R. But how can I calculate inv(H) if each element of H is a matrix? I cannot simply concatenate these smaller matrices into a bigger matrix H and then invert it - since this will give a different result than, e.g. inverting matrix H with symbolic values and then substituting these symbolic values with smaller matrices (because of the pointwise product present in the set of equations).

Up to now I came up with 1 solution. I will create the matrix H with 20x20 symbolic values, I will invert it, and then I will evaluate each cell of the resulting inverted matrix with 'subs'. 

H = sym('A',[20 20]);
invmat = inv(H) ;
% here I load 400 smaller matrices with appropriate names
invmat_11 = subs(invmat(1,1));

But inversion of such a matrix is to complicated to be calculated on any medium class computer so I never managed to run this code. Do you know any other way of calculating an inversion of matrix H - or directly solving for S?

I was asked for some simple example: Consider equation R = H*S (S is unknown). Let's say H=[A B; C D], where A,B,C,D are 2x2 matrices, e.g. A = [A11 A12; A21 A22]

and R and S are 2x1 matrices, e.g.

R = [R1;R2]

To calculate for S I need to solve inv(H)*R, symbolically inv(H) =

[ -D/(B*C - A*D), B/(B*C - A*D)] [ C/(B*C - A*D), -A/(B*C - A*D)]

Now I can substitute A,B,C and D with real 2x2 matrices and calculate the inversion of H:

inv(H) = [H1 H2; H3 H4] where

H1 = -D/(B*C - A*D)

This constitutes calculation of inv(H). Now I need to multiply inv(H) with R (to solve for S):

S1 = H1*R1 + H2*R2 S2 = H3*R1 + H4*R2

but please note, that all H1 to H4 and R1 to R2 are matrices, and * symbolizes pointwise product.

pewter_cauldron
  • 73
  • 6
  • 2
    Could you add your (slow) solution to the question. Maybe the performance can be increased biased on that solution. How large are your smaller matrices? – Daniel Jul 23 '15 at 09:05
  • 1
    Just a note - it is not efficient to calculate `inv(H)*R` (and even fails if H is rank deficit). The usual solution is to solve the equation - eg. use [mldivide](http://de.mathworks.com/help/symbolic/mldivide.html) - which is short written as `H\R` – bdecaf Jul 23 '15 at 09:10
  • @Daniel, my slow solution is already in the question: my computer never managed to compute H = sym('A',[20 20]); invmat = inv(H); I even tried to apply the blockwise inversion ( https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion ) so that I would need to calculate only 10x10 matrix inversion - but that is still to memory-consuming. – pewter_cauldron Jul 23 '15 at 09:34
  • @bdecaf I know that mldivide is much faster, but I don't now how to apply it when pointwise product is present in the equations. It would work if I concatenate all smaller matrices into a bigger matrix H but then the results would be wrong. – pewter_cauldron Jul 23 '15 at 09:37
  • 1
    Could you add a simple example (say 2x2 matrices) with expected output? – bdecaf Jul 23 '15 at 09:42
  • 1
    @bdecaf I added some example. – pewter_cauldron Jul 23 '15 at 10:22
  • So just for completeness R1, A1 and S1 are 2x2 matrices? Would the answer be correct if you replace in the H1 and S1 formula the multiplication and division with the pointwise variants? – bdecaf Jul 23 '15 at 12:36
  • @bdecaf I believe so. If, when calling Matlab inv function, Matlab would pointwise multiply and pointwise divide the elements, I think the results would be ok. But the problem will still be extremely compunational-consuming. – pewter_cauldron Jul 23 '15 at 13:38
  • so `ldivide` and `rdivide` would be ok? – bdecaf Jul 23 '15 at 13:46
  • Well, technically yes, but it would have to be applied to blocks of matrix H, which I don't think that is possible. And again: my medium class computer is unable to solve such a problem so I'm also looking for some other way of solving this. – pewter_cauldron Jul 23 '15 at 13:56
  • Do I understand correctly that you actually don't need the symbolic result, but use it as intermediate step? I believe you could write it down in tensor form. Then you could use the `kron` command to get the corresponding matrix which can be solved numerically. – bdecaf Jul 26 '15 at 10:50

1 Answers1

1

I found the best solution to solve this set of equations. Actually, it is embarrassingly easy: one just has to notice that these equations can be rewritten in the form:

first_element_of_r1 = first_element_ofs1*first_element_of_h1_1 + ...

and this is due to presence of piecewise product in the equations. Now each element of r1 to r20 matrices can be solved independently in a loop (or parallel loop). Thanks everyone for trying to help me.

pewter_cauldron
  • 73
  • 6