this is purely out of interest; how would one go about solving Ax = b (A nxn matrix, x,b vectors) for A, knowing that x,b are such that A is at least symmetric positive definite.

I guess for starters could restrict it to A being 2x2, still not sure how to go about it.

By hand you'd just use Gaussian elimination or Jordan elimination. If you aren't familiar with it, go to the library or your favorite math teacher and ask them for a book on linear algebra.


With only that brief description, it looks an awful lot like homework from a 2nd-year linear algebra class. When a person is able to throw around words like a "symmetric positive definite" matrix, either they've got enough math background to understand it or they're working from a college linear algebra book.

With only that brief description, it looks an awful lot like homework from a 2nd-year linear algebra class.

That was my first impulse too -- I was going to point him to the Cholesky factorization (which is the fastest way I'm aware of to invert positive definite matrices) -- but then I noticed that he wasn't trying to solve for x, but rather to find some A that would solve his underdetermined problem.

this is purely out of interest; how would one go about solving Ax = b (A nxn matrix, x,b vectors) for A, knowing that x,b are such that A is at least symmetric positive definite.

I guess for starters could restrict it to A being 2x2, still not sure how to go about it.

There isn't going to be a unique answer. You have a linear constraint on A, and a convex constraint (that it be positive definite), so you've restricted it to some convex set, but that's it; you haven't specified it uniquely.

The easiest way to find such an A, if you have the right optimization libraries, is to treat this as a semidefinite program. You could solve the feasibility problem,

min[sub]A,z[/sub] z
s.t.
A + z I > 0
A[sup]T[/sup] = A
A x = b .

If z[sub]opt[/sub] in your optimal pair (A[sub]opt[/sub],z[sub]opt[/sub]) is negative, your problem is feasible, and A[sub]opt[/sub] is an answer.

Typically you'd solve such a problem with an interior point method. I usually use cvx for this, which is a MATLAB frontend for the MATLAB library SeDuMi. I'm not sure what the best pure-C/C++ library is for this, but you can find a list of solvers at the Wikipedia article.

thanks Emergent, i had a feeling that it would not in general be unique; I don't actually have a situation in which I need to do this, like I said it was purely out of the interest but it's interesting to see that solving it is not a trivial thing to do; it seems like there 'should' be some algebraic trick to simplify the problem but one might have said the same thing w.r.t. some of the incredibly difficult to prove but simple looking problems.

Here's another way to think about this:

Since A is positive definite and symmetric, by the spectral theorem it factors as A = R^T D R, where R is a rotation matrix and D is a positive diagonal matrix. Consequently, the equation

A x = b

can be rewritten

D R x = R b .

Then, once you have R, it's easy to find D, and it will be positive so long as R x and R b have the same signs elementwise. Geometrically, this means R x and R b are in the same orthant, and without loss of generality it might as well be the positive orthant. So we've reduced the problem to finding a rotation that brings two given vectors into the positive orthant.

This is not always possible, and the easiest example is b = -x. A matrix 'A' fails to be positive definite if there is a vector 'z' s.t. z^T A z <= 0. Letting z = x, you get z^T A z = -||x||^2 < 0. Sometimes, your equation has no solution.

In general, you can find a rotation that does this only if the angle between the two vectors is not too large. Specifically, they'll fit into the positive orthant IFF the angle between them is less than 90 degrees. In other words, their inner product needs to be positive;

<x, b> > 0.

This is nice, because it characterizes in simple terms when your equation has a solution.

Next, to build 'A,' we need to actually construct a rotation that takes x and b to the positive orthant. Letting

v1 = (1,0,...,0)
v2 = (0, 1/sqrt(n), ..., 1/sqrt(n))

m = normalize(x/||x|| + b/||b||)
d = normalize(x/||x|| - b/||b||)

vm = normalize (v1 + v2)
vd = normalize(v1 - v2)

I would choose a rotation that maps,

m -----> vm
d -----> vd.

To do this, build
1.) a matrix R1 whose first two columns are m and d, and whose subsequent columns are random and orthonormal to all previous ones (do this in a Gram-Schmidt-like-way)
2.) a matrix R2, in the same style as R1, whose first two columns are vm and vd;
then let R = R2 R1^T.

Now that you have R, the diagonal matrix D is easy to find, and then you can construct A as A = R^T D R.

There you go! Conditions on when a solution exists, and constructive way to produce one when it does.