Solving an Incomplete Matrix such that A*transpose(A) = eye(n)

Hello,

I'd like to write a function that takes in a matrix composed of numbers and variables, and solves for the value of the variables given the relationship A * transpose(A) = eye(n) where eye(n) is an n*n identitiy matrix. I've gotten to the point where I can get a system of equations in matrix form, but I'm having trouble extracting the variables from that system. Any and all help would be greatly appreciated!

syms a b c;
A = [0.2193 0.1864 0.9589; a (-0.0160) (-0.2147);b c (-0.1845)];
B = transpose(A);
C = (A*B);
I = eye(3);
A = C == I;

% Given this relationship, I'd now like to solve for a, b, and c. How can I do this?

Answer 1

[John D'Errico]

There exist NO such values for a,b,c such that A*A' == eye(3).

Look carefully at A*A'.

vpa(A*A')
ans =
[                       1.00232666,      0.2193*conj(a) - 0.20885823, 0.2193*conj(b) + 0.1864*conj(c) - 0.17691705]
[            0.2193*a - 0.20885823,           a*conj(a) + 0.04635209,       a*conj(b) - 0.016*conj(c) + 0.03961215]
[ 0.2193*b + 0.1864*c - 0.17691705, b*conj(a) - 0.016*c + 0.03961215,           b*conj(b) + c*conj(c) + 0.03404025]

What is the (1,1) element of that matrix? 1.0023266. Which must be exactly 1, yet that is not a function of a, b or c at all.

Next, consider the (2,1) element. Here we will learn that

0.2193*a - 0.20885823 == 0

which implies

a = 0.20885823/0.2193 = 0.95239

But then we see from the (2,2) element that

a*a' + 0.04635209
ans =
      0.95339

which must be 1, not 0.95...

I can go on, but as I said, there are no values for a, b, and c that will yield the identity matrix, at all closely. There are cetainly no such values that will be accurate to 4 significant digits.

One important thing to learn from this is you should strongly avoid those short, 4 digit approximations to numbers. If you were to perturb the numeric values in A, you might be able to find a solution.

An alternative is to identify IF there is some set of valiues for a,b,c such that the product is as close as possible to eye(3).

Afun = @(a,b,c) [0.2193 0.1864 0.9589; a  -0.0160  -0.2147 ;b c -0.1845];
obj = @(abc) Afun(abc(1),abc(2),abc(3))*Afun(abc(1),abc(2),abc(3))' - eye(3);
abc = lsqnonlin(obj,[1 1 1])

abc = lsqnonlin(obj,[1 1 1])

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>
abc =
      0.97595     -0.02467      0.98244

So a solution was found. How close did we come to a valid solution?

obj(abc)
ans =
    0.0023267     0.005168   0.00080037
     0.005168   -0.0011659  -0.00018376
   0.00080037  -0.00018376  -0.00015485

If that was all zeros, to within +/-0.0001, then I would say a solution reasonably exists, and the problem was merely that you gave only approximate values for the constant terms in your matrix. Since those errors are too large by as much as a factor of 50:1, that suggests those numbers are inaccurate by more than just the rounding error out at the 4th decimal digit.

I suppose, with some additional effort, we could next find the smallest perturbations to the original matrix A such that an exact solution exists.

[Matthew Kolodzik]

Thank you! I see what you're saying now is true.

[Matthew Kolodzik]

Just out of curiousity, why is it that you can't generalize this to be solved in matlab, but it is posssible to break this down into a system of equations by hand and solve for the values of a, b, and c in that manner?

[John D'Errico]

That is a valid question. When you formulate a matrix equality like that, there are 9 elements of the equality. They must all be exactly satisfied for a tool like solve to work. So we have what is effectively 9 equations. But your constants are not exact. So what happens is you will never find 3 parameters that will solve all 9 equations EXACTLY. MATLAB's solve cannot handle that, because solve is an exact solver. It does not understand least squares things, thus finding the solution that minimizes the errors of all 9 elements in an overall sense.

However, suppose we pick 3 of the equations. ANY 3 will probably suffice, as long as they are independent, and are sufficient to solve for each of the unknowns. That is, if one of the equations reduces to a constant == another constant, then it is useless. Or, if all three elements reduce to terms that have only a and b in them, but not c, then again you are lost.

But as long as we pick 3 of those 9 possible choices in an intelligent way, then we can compute a, b, and c. The problem arises if we make a different choice of which equations? Remember, these "known" constants were not exactly known. You only gave them to us with 4 significant digits. So we will get slightly different values, depending on how we make our choices. Here is the matrix you gave us.

syms a b c
A = [0.2193 0.1864 0.9589; a  -0.0160  -0.2147 ;b c -0.1845]
A =
[ 2193/10000, 233/1250,  9589/10000]
[          a,   -2/125, -2147/10000]
[          b,        c,   -369/2000]

Now compute the matrix problem.

M = vpa(A*A' - eye(3))
M =
[                       0.00232666,      0.2193*conj(a) - 0.20885823, 0.2193*conj(b) + 0.1864*conj(c) - 0.17691705]
[            0.2193*a - 0.20885823,           a*conj(a) - 0.95364791,       a*conj(b) - 0.016*conj(c) + 0.03961215]
[ 0.2193*b + 0.1864*c - 0.17691705, b*conj(a) - 0.016*c + 0.03961215,           b*conj(b) + c*conj(c) - 0.96595975]

Our goal, IF this is satisfied, is to have ALL 9 elements of that matrix be zero. The first one is like trying to get blood from a rock. As hard as I try, nothing will suffice.

M(1,1) 
ans =
0.00232666

I can try setting M(2,1) to zero, as I could have done with M(1,2). They are the same, since either will yield the same solution for a.

M(2,1)
ans =
0.2193*a - 0.20885823

M(1,2)
ans =
0.2193*conj(a) - 0.20885823

solve(M(1,2) == 0,a)
ans =
0.9523859097127222982216142

But would M(2,2) have given us the same value for a?

M(2,2) == 0
ans =
a*conj(a) - 0.95364791 == 0

solve(M(2,2) == 0,a)
ans =
 -0.9765489798264089119062242
  0.9765489798264089119062242

In fact, that yields two solutions, neither of which was the same as if we had started in a different place.

Regardless of how we choose to solve for a here though, we could now pick another element of M. Perhaps

M(3,2) == 0
ans =
b*conj(a) - 0.016*c + 0.03961215 == 0

Since we know a now we can find c. And then we could recover b from:

M(3,1) == 0
ans =
0.2193*b + 0.1864*c - 0.17691705 == 0

But you do need to recognize that if I chose things differently, I would have gained a subtly different solution, by choosing a different route.

The symbolic toolbox cannot resolve all of the various possible solutions I might find using the one at a time approach. I was able to find an overall solution using lsqnonlin, because lsqnonlin looks for a least squares solution. It tries to minimize ALL of the errors of each of those 9 elements at once. And that may more reasonably have an overall solution that minimizes the total error. So effectively, I was able to solve it in MATLAB, as long as I chose to solve it the right way.

Answer 2

[Sai Sri Pathuri]

You may add the below line

s = solve(A,[a,b,c]);

The solution can be obtained as

aValue = s.a;
bValue = s.b;
cValue = s.c;

The above command returns empty structure when there is no solution

[John D'Errico]

While the call to solve would certainly be valid under some circumstances...

You should recognize the solve as you have suggested creates effectively a system of 9 equations, in 3 unknowns. As such, solve on the overdetermined problem will ALWAYS fail to find a solution. Worse, since one of those equations is the constant one, where we need to solve

1.0023266 == 1

the solution must always fail, even if it was not overdetermined?

If you suggest using solve, you should point out that it must fail, not that it MAY fail, by the hint that it returns an empty result if there was no solution. This is because someone would likely want to then use solve, and they will certainly be surprised at the result.

[Matthew Kolodzik]

I tried this but what you said about there being no solution was what was returned.

[Matthew Kolodzik]

John,

Thank you so much for that explaination. It's very thorough and answeres all my questions. Thanks so much!