0

Given a quadratic form $q: E \to \mathbb{R}$ (where $E$ is some vector space), it's always possible to find a decomposition of the form $$q(v) = \sum_{i=1}^n \alpha_i l_i(v)^2,$$ where the $\alpha_i \in \mathbb{R}$ are scalars and $l_i : E \to \mathbb{R}$ are independent linear forms. In French it's called the Gaussian decomposition ("décomposition de Gauss") of $q$; I don't know if this has a standard name in English.

Is there a way to find these decompositions easily in Mathematica? One way would be, given a symmetric square matrix $M$, to find an invertible matrix $P$ such that $P^T M P$ is diagonal; the matrix $P$ would represent the predual basis of $(l_i)$). As far as I can tell, Mathematica only has Orthogonalize, which only works for positive definite quadratic forms.

To give an example of what I'm asking. Say that $q : \mathbb{R}^3 \to \mathbb{R}$ is defined as $$q(x,y,z) = 2 x^2 + 12 x y + 17 y^2 + 4 x z + 16 y z - z^2.$$ Then, completing the squares, one gets $$q(x,y,z) = 2 (x + 3 y + z)^2 - (y - 2 z)^2 + z^2,$$ so I would either want an output of the form {{2,-1,1}, {{1,3,1}, {0,1,-2}, {0,0,1}}} (for the dual basis) or {{2,-1,1}, {{1, -3, -7}, {0, 1, 2}, {0, 0, 1}}} (for the predual basis).

Note that Orthogonalize fails, as q is not positive definite:

q[{x_, y_, z_}] = 2  (x + 3 y + z)^2 - (y - 2 z)^2 + z^2
b[u_, v_] := (q[u + v] - q[u] - q[v])/2
Orthogonalize[IdentityMatrix[3], b]
(* Orthogonalize::ornfa: The second argument b is not an inner product function, which always should return a number or symbol. If its two numeric arguments are the same, it should return a non-negative real number. *)
Najib Idrissi
  • 576
  • 3
  • 11

1 Answers1

2

"Eigensystem" is the function you are looking for. Here is an example:

Given some symmetric square matrix:

n = 3;
m = RandomReal[{-1, 1}, {n, n}]; m += Transpose[m];

The eigenvalues (es[[1]]) and eigenvectors (es[[2]]) are:

es = Eigensystem[m];

Now the different linear forms can be written:

lf[x_, i_] := es[[1, i]]  (es[[2, i]] . x)^2

And with this we get the quadratic form:

qf[x_] = Sum[lf[x, i], {i, n}];

Now let us test this:

t = RandomReal[{-1, 1}, n];
t . m . t == qf[t]

True
Daniel Huber
  • 51,463
  • 1
  • 23
  • 57
  • Very strange. If I test your code on random real matrices, it does appear to work. But if I test on my example (m = {{2, 6, 2}, {6, 17, 8}, {2, 8, -1}}) I get wrong results. And if I replace RandomReal by RandomInteger (say between -10 and 10), your code also doesn't work. Besides, I'm not able to see why your answer is sound mathematically. Could you explain how you came up with it? – Najib Idrissi Feb 29 '24 at 12:27
  • As far as I can tell, what you're doing is the following. Supposed $Mv_i = \lambda_i v_i$ is the eigensystem, and $(v_i^)$ is a dual basis for the usual inner product (i.e. $v_i \cdot v_j^ = \delta_{ij}$). Then for all vectors, $u = \sum_i (v_i \cdot u) v_i^$, and so $q(u) = \sum_i (v_i \cdot u)^2 q(v_i^)$. The map $qf$ you're defining is rather $\sum_i (v_i \cdot u)^2 \lambda_i$, so, essentially, you seem to be saying that $\lambda_i = q(v_i^) = v_i^ \cdot (Mv_i^*)$. Why is this true? – Najib Idrissi Feb 29 '24 at 12:41
  • Based on some experimentation, it seems that Eigensystem uses a different algorithm when the matrix is exact or approximate. (When m is exact, all the eigenvectors have their last coordinate equal to 1, for starters.) My guess is that when the matrix is approximate, then the eigenvectors it finds have a certain structure that makes the above work, but when it's exact, it breaks down. – Najib Idrissi Feb 29 '24 at 12:44
  • 1
    The mystery comes from the norm of the eigenvectors. My methods needs normalized eigenvectors. Eigensystem normalizes the eigenvectors for real matrices. However, for accurate numbers, the eigenvectors are NOT normalized. If you insert es[[2]] = Normalize /@ es[[2]];to normalize the vectors, the code will also work with integer matrices. – Daniel Huber Feb 29 '24 at 13:14