Solve the vector-matrix equation. Minimize the length of the desired n-dimensional vector

Question

There is the following vector-matrix equation:

$$\mathbf x^\top\mathbf M\mathbf x=\begin{bmatrix}x_1&x_2&x_3\end{bmatrix}\begin{bmatrix}a&b&c\\d&e&f\\g&h&i\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=C$$

where $C$ is a nonzero constant, $\begin{bmatrix}x_1&x_2&\cdots&x_n\end{bmatrix}^\top$ is a vector of real numbers, and $\mathbf M$ is a matrix of real numbers of dimension $n \times n$.

It is required to find a vector $\mathbf x$ that has a minimum Euclidean norm, i.e.

$$\min_{\mathbf x^\top\mathbf M\mathbf x=C}\|\mathbf x\|_2$$

1. How to solve this problem analytically within the framework of the vector-matrix apparatus?
2. How to solve this problem with optimization in Mathematica?

Answer 1

Really this should have a concrete example in Mathematica copy-pastable format. Anyway, I'll show a few ways using a made-up example. We'll eventually get it down to straight linear algebra.

To keep it simple I'll assume the matrix has distinct eigenvalues. This restriction can be lifted of course. Also as has been noted in comments, without loss of generality we can assume the matrix is symmetric, since adding the transposed equation gives us an equivalent problem with a symmetric matrix. It should be mentioned that the problem need not have a solution, e.g. if the matrix is positive definite and the constant is negative. Also, multiplying through by -1 if need be, we can assume the constant is positive. Finally note that the solution is only unique up to sign.

We'll start with an example.

SeedRandom[1234];
c = 3;
vars = {x, y, z};
obj = vars.vars;
mat = RandomInteger[{-10, 10}, {3, 3}];
symmat = mat + Transpose[mat]

(* Out[28]= {{-18, 1, 17}, {1, -20, 0}, {17, 0, -10}} *)

The most pedestrian thing we can do is to just invoke a minimizer directly on objective and constraint.

NMinimize[{obj, vars.symmat.vars == c}, vars]

(* Out[29]= {0.861901, {x -> -0.576638, y -> -0.024558, z -> -0.727177}} *)

We can say more than this however. In 3D (actually in general) what we are looking for is the smallest radius sphere centered at the origin that intersects with the quadric surface defined by a given level set from the quadratic form defined by the symmetric matrix. A standard thing to do is use a Lagrange multiplier formulation.

solns = 
  NSolve[
    Flatten[{vars.symmat.vars == c, 
      Thread[Grad[obj, vars] == 
        lambda*Grad[vars.symmat.vars, vars]]}], 
    Join[{x, y, z}, {lambda}]];
realsolns = Select[vars /. solns, FreeQ[#, Complex] &]

(* Out[53]= {{-0.576638, -0.024558, -0.727177}, {0.576638, 0.024558, 
  0.727177}} *)

So same solution(s).

As has also been noted, one can diagonalize the system. The diagonalizing matrix will have orthogonal components hence can be made orthogonal my normalizing rows if necessary (for approximate valued matrices, Eigensystem does this automatically). Once we do this things get simpler. The objective function does not change since the Euclidean metric does not change under orthogonal transformation.

We diagonalize, and solve in the rotated system.

{evals, evecs} = Eigensystem[N@symmat];
{min, rvals} = 
 NMinimize[{obj, {evals.vars^2 == c, Thread[{x, y, z} >= 0]}}, {x, y, 
   z}]

(* Out[57]= {0.861901, {x -> 0., y -> 0., z -> 0.928386}} *)

Rotate back to recover the solution in the original coordinate system.

vals = Transpose[evecs].(vars /. rvals)

(* Out[58]= {0.576638, 0.024558, 0.727177} *)

Again the same solution (good). But now we see that the diagonalized system has a "nice" solution: it lies on a coordinate axis. That's because the primary axes of the diagonalized quadratic form are just the standard Cartesian axes. This makes things easy. We can, for example, replace the squares of the variables by new variables, reducing to linear programming (the new variables have to be positive so we'll get real square roots on back-transforming). We show this by recovering the intermediate solution above (the one prior to rotating back to the original system).

{min, svals} = 
 NMinimize[{Total[vars], {evals.vars == c, Thread[vars >= 0]}}, vars]
Sqrt[vars /. svals]

(*Out[62]= {0.861901, {x -> 0., y -> 0., z -> 0.861901}}

Out[63]= {0., 0., 0.928386} *)

But really this is all overkill. We should just observe that the minimizer has to be on a coordinate axis, and it will be the axis with the largest positive eigenvalue: coordinates corresponding to negative eigenvalues just add to the magnitude, and the ellipsoidal part of the quadric surface level set in the subspace corresponding to positive eigenvalues will be thinnest along the axis of the largest one. So we should just look at the eigenvalues.

evals

(* Out[64]= {-31.5177, -19.963, 3.48068} *)

Our only candidate is the third one since the others are negative. A bit of reasoning will show that the value for the corresponding coordinate is as below.

Sqrt[c/evals[[3]]]

(* Out[68]= 0.928386 *)

So there it is. The problem can be solved by (1) Diagonalizing the quadratic form. (2) Setting all but one coordinate to zero and the coordinate corresponding to the largest eigenvalue to the square root of the constant divided by that eigenvalue. (3) Use the conversion matrix from the diagonalization to rotate back, that is, get a solution in the original system. The other solution comes from negating the solution in step (2). All this basically just elaborates on the method posted in a comment by @mikado.

There is another thing worthy of note. The Lagrange multipliers in the system of equations are just reciprocals to the eigenvalues of the rotated system.

Sort[Union[lambda /. solns, SameTest -> (Norm[#1 - #2] < 10^(-4) &)]]

(* Out[82]= {-0.0500926, -0.0317282, 0.2873} *)

Sort[1/evals]

(* Out[83]= {-0.0500926, -0.0317282, 0.2873} *)

I sort of recall first encountering this connection of minimal length solutions with quadratic form constraints, and Lagrange multipliers, end eigenvalues, way back in grad school.