Solve matrix equation in matlab

linear-algebra, matlab, matrix

Solution

Let `X` be the collection of all 100,000 `x`s you got (such that the `i`-th column of `X` equals the `x_i`-th vector). In the same manner we can define `Y` and `C` as 2D collections of `y`s and `c`s respectively.

What you wish to solve is for `A` and `B` such that

C = AX + BY

You have 2 * 50,000^2 unknowns (all entries of `A` and `B`) and `numel(C)` equations.

So, if the number of data vectors you have is 100,000 you have a single solution (up to linearly dependent samples). If you have more than 100,000 samples you may seek for a least-squares solution.

Re-writing:

C = [A B] * [X ; Y]  ==>  [X' Y'] * [A';B'] = C'

So, I suppose

[A' ; B'] = pinv( [X' Y'] ) * C'

In matlab:

ABt = pinv( [X' Y'] ) * C';
A = ABt(1:50000,:)';
B = ABt(50001:end,:)';

Correct me if I'm wrong...

EDIT: It seems like there is quite a fuss around dimensionality here. So, I'll try and make it as clear as possible.

Model: There are two (unknown) matrices `A` and `B`, each of size 50,000x50,000 (total 5e9 unknowns). An observation is a triplet of vectors: (`x`,`y`,`c`) each such vector has 50,000 elements (total of 150,000 observed points at each sample). The underlying model assumption is that an observation is generated by `c = Ax + By` in this model. The task: given `n` observations (that is `n` triplets of vectors { (`x_i`, `y_i`, `c_i`) }_`i=1..n`) the task is to uncover `A` and `B`.

Now, each sample (`x_i`,`y_i`,`c_i`) induces 50,000 equations of the form `c_i = Ax_i + By_i` in the unknown `A` and `B`. If the number of samples `n` is greater than 100,000, then there are more than 50,000 * 100,000 ( > 5e9 ) equations and the system is over constraint.

To write the system in a matrix form I proposed to stack all observations into matrices:

- A matrix `X` of size 50,000 x `n` with its `i`-th column equals to observed `x_i`

- A matrix `Y` of size 50,000 x `n` with its `i`-th column equals to observed `y_i`

- A matrix `C` of size 50,000 x `n` with its `i`-th column equals to observed `c_i`

With these matrices we can write the model as:

C = A*X + B*Y

I hope this clears things up a bit.

Thank you @Dan and @woodchips for your interest and enlightening comments.

EDIT (2): Submitting the following code to octave. In this example instead of 50,000 dimension I work with only 2, instead of `n=100,000` observations I settled for `n=100`:

n = 100;
A = rand(2,2);
B = rand(2,2);
X = rand(2,n);
Y = rand(2,n);
C = A*X + B*Y + .001*randn(size(X)); % adding noise to observations 
ABt = pinv( [ X' Y'] ) * C';

Checking the difference between ground truth model (`A` and `B`) and recovered `ABt`:

ABt - [A' ; B']

Yields

  ans =

   5.8457e-05   3.0483e-04
   1.1023e-04   6.1842e-05
  -1.2277e-04  -3.2866e-04
  -3.1930e-05  -5.2149e-05

Which is close enough to zero. (remember, the observations were noisy and solution is a least-square one).

Problem

I have an equation of the type `c = Ax + By` where `c`, `x` and `y` are vectors of dimensions say 50,000 X 1, and `A` and `B` are matrices with dimensions 50,000 X 50,000. Is there any way in Matlab to find matrices `A` and `B` when `c`, `x` and `y` are known? I have about 100,000 samples of `c`, `x`, and `y`. `A` and `B` remain the same for all.

Original source