2
$\begingroup$

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?
$\endgroup$
  • 2
    $\begingroup$ Note that you can replace M with a symmetrical matrix that gives the same result. Your solution is given by the eigenvector corresponding to the largest eigenvalue. You might need to think further if C could be negative. $\endgroup$ – mikado Oct 12 at 12:28
  • $\begingroup$ Hello, mikado. Thank you for your response. Did i understand you correctly that i need to find the eigenvectors of the matrix M and choose from them the one that corresponds to the largest eigenvalue? I would like to understand why this is so, because am i new to this section of mathematics? Will this be a global solution? P.S. Yes, C can be negative - why should this be taken into account? $\endgroup$ – Andy Sol Oct 12 at 13:17
  • $\begingroup$ This can be formulated as a a minimization problem with quadratic objective function and quadratic equality constraint(s). I do not think that this can be boiled down to a single linear system of equations. But it can be solved, e.g. by Newton's methos. You might want to use FindMinimum for that. $\endgroup$ – Henrik Schumacher Oct 12 at 14:07
  • $\begingroup$ I have already solved this problem with the help of Mathematica optimization modules. But, i want to get an analytical solution. $\endgroup$ – Andy Sol Oct 12 at 14:11
  • 2
    $\begingroup$ In outline, a symmetric matrix can be diagonalised with an orthogonal matrix. You can easily solve the problem for a diagonal matrix. You then rotate the solution with the orthogonal matrix (giving the eigenvector mentioned above). $\endgroup$ – mikado Oct 12 at 14:35
9
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ It can be used for varying C[t]. Best I think would be to find where the eigenvalues become multiple since those values of t will demarcate which eigenvalue/position becomes largest (they can only cross at points of multiplicity). I guess M[t] can also vary for that matter. Might make the equation solving more complicated though. An alternative is to just compute for specific values of t, so the matrix and value are effectively constants. $\endgroup$ – Daniel Lichtblau Oct 12 at 20:16
  • 1
    $\begingroup$ (1) Depends on exactly what you want, but certainly most can be automated. I was a bit thick about what has to be done though. For the constant it is just a matter of checking values of t for which it changes sign, since that determines the eigenvalues that are important. $\endgroup$ – Daniel Lichtblau Oct 13 at 14:15
  • 1
    $\begingroup$ (2) For eigenvalue multiplicity, it is the discriminant of the characteristic polynomial, as a function of the eigenvalue variable (call it lam) that matters. Programatically, that's discrim=Discriminant[CharacteristicPolynomial[M[t],lam],lam]. This gives a function of t, the zeros of which give points of multiplicity for the eigenvalues. So you'd find those values of t via Solve[discrim==0,t]. $\endgroup$ – Daniel Lichtblau Oct 13 at 14:15
  • 1
    $\begingroup$ (3) Last, you'll want to know when eigenvalues are zero since their sign matters as well. That's simply Solve[Det[M[t]]==0,t]. Now you have a bunch of important values of t, where all possible trnsitions can occur in terms of which eigenvalue is the largest of same magnitude as C[t]. $\endgroup$ – Daniel Lichtblau Oct 13 at 14:15
  • 1
    $\begingroup$ I should add that this all really belongs in the question, along with a smallish representative example. Trying to work out the details in comments is both tricky and to some extent in violation of the charter of the forum (as the moderators will soon enough tell us, I suspect). $\endgroup$ – Daniel Lichtblau Oct 13 at 14:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.