4
$\begingroup$

For given $n\times n$ real symmetric positive definite matrices $X_{1}, X_{2}$, let $Y := X_{1}^{-1}X_{2}$, and let $I$ be the identity matrix.

I would like to solve the following equation for the unknown $p>0$, in terms of functions of the matrices $X_{1}$ and $X_{2}$:

$$\text{trace}\left[\left(I + pY\right)^{-1} \left(I - p^{2}Y\right)\right] = 0.$$

I tried expanding the inverse in Neumann series but I am not sure how that helps to solve for $p$. This is not homework.

$\endgroup$
10
$\begingroup$

Write $Y=SDS^{-1}$, where $D$ is diagonal (since $X_1$ and $X_2$ are psd, $Y$ is diagonalizable). Then, observe that

\begin{equation*} f(p) = \text{tr}(S^{-1}(I+pSDS^{-1})^{-1}SS^{-1}(I-p^2SDS^{-1})S). \end{equation*} This simplifies into a simpler function of $p$ involving only the diagonal values in $D$. That particular equation you can solve using Newton's method.


EDIT: as noted in my comments below, the same idea also works more generally if we use the Schur decomposition to reduce matrices to an upper triangular form; we do not need to assume diagonalizability in this case.

$\endgroup$
4
  • $\begingroup$ Why is the first statement true (that $Y$ is diagonalizable)? $\endgroup$ – Igor Rivin Jan 18 '17 at 4:33
  • 3
    $\begingroup$ @IgorRivin: Since $X_{1}$ is symmetric positive definite, so is $X_{1}^{-1}$. Then we have this $\endgroup$ – Abhishek Halder Jan 18 '17 at 5:50
  • $\begingroup$ @AbhishekHalder Cool, I did NOT know this trick! $\endgroup$ – Igor Rivin Jan 18 '17 at 14:21
  • 1
    $\begingroup$ using Rodrigo's observation and my comment to it, the above idea extends to general matrices by reducing them to upper triangular using Schur decomposition, and then solving an equation involving eigenvalues of the matrix $Y$ (which need not be diagonalizable any more). $\endgroup$ – Suvrit Jan 18 '17 at 22:51
4
$\begingroup$

Given $\mathrm A \in \mathbb R^{n \times n}$, we define $f : \mathbb R \to \mathbb R$ as follows

$$f (x) := \mbox{tr} \left( (\mathrm I_n + x \mathrm A)^{-1} (\mathrm I_n - x^2 \mathrm A) \right)$$

It is not necessary that $\mathrm A$ be diagonalizable. Using the Schur decomposition $\mathrm A = \mathrm Q \mathrm U \mathrm Q^{\top}$, where $\mathrm Q$ is orthogonal and $\mathrm U$ is upper triangular with the eigenvalues of $\mathrm A$ on its main diagonal, we have

$$\begin{array}{rl} f (x) &= \mbox{tr} \left( (\mathrm I_n + x \mathrm A)^{-1} (\mathrm I_n - x^2 \mathrm A) \right)\\ &= \mbox{tr} \left( (\mathrm Q \mathrm Q^{\top} + x \, \mathrm Q \mathrm U \mathrm Q^{\top})^{-1} (\mathrm Q \mathrm Q^{\top} - x^2 \, \mathrm Q \mathrm U \mathrm Q^{\top}) \right)\\ &= \mbox{tr} \left( \mathrm Q (\mathrm I_n + x \, \mathrm U )^{-1} \mathrm Q^{\top} \mathrm Q (\mathrm I_n - x^2 \, \mathrm U) \mathrm Q^{\top} \right)\\ &= \mbox{tr} \left( (\mathrm I_n + x \, \mathrm U)^{-1}(\mathrm I_n - x^2 \, \mathrm U) \right)\end{array}$$

where both $\mathrm I_n + x \, \mathrm U$ and $\mathrm I_n - x^2 \, \mathrm U$ are upper triangular. Note that

  • The inverse of an upper triangular matrix is also upper triangular. Since $\mathrm I_n + x \, \mathrm U$ is upper triangular, $(\mathrm I_n + x \, \mathrm U)^{-1}$ is also upper triangular.

  • The product of two upper triangular matrices is also upper triangular with the main diagonal being the entrywise product of the main diagonals of the factors. Hence, we can conclude that $(\mathrm I_n + x \, \mathrm U)^{-1}(\mathrm I_n - x^2 \, \mathrm U)$ is also upper triangular.

Suppose that $\mathrm A$ has $m$ distinct eigenvalues with multiplicities $n_1, n_2, \dots, n_m$. Hence,

$$f (x) = \sum_{k=1}^m n_k \left( \dfrac{1 - x^2 \, \lambda_k (\mathrm A)}{1 + x \, \lambda_k (\mathrm A)} \right)$$

Thus, the equation $f (x) = 0$ can be written as follows

$$n_1 \left( \dfrac{1 - x^2 \, \lambda_1 (\mathrm A)}{1 + x \, \lambda_1 (\mathrm A)} \right) + n_2 \left( \dfrac{1 - x^2 \, \lambda_2 (\mathrm A)}{1 + x \, \lambda_2 (\mathrm A)} \right) + \cdots + n_m \left( \dfrac{1 - x^2 \, \lambda_m (\mathrm A)}{1 + x \, \lambda_m (\mathrm A)} \right) = 0$$


Example

Suppose we have

$$\mathrm A = \begin{bmatrix} 1 & 0 & 0\\ 0 & 2 & 1\\ 0 & 0 & 2\end{bmatrix}$$

which is in Jordan normal form and, thus, is already upper triangular. Matrix $\mathrm A$ has $m = 2$ distinct eigenvalues with multiplicities $n_1 = 1$ and $n_2 = 2$. Using SymPy:

>>> from sympy import *
>>> A = Matrix([[1,0,0],[0,2,1],[0,0,2]])
>>> x = Symbol('x')
>>> ((eye(3) + x * A)**-1 * (eye(3) - x**2 * A)).trace()
2*(-2*x**2 + 1)/(2*x + 1) + (-x**2 + 1)/(x + 1)

We obtain the equation

$$\frac{1 - x^{2}}{x + 1} + 2 \left(\frac{1 - 2 x^{2}}{2 x + 1}\right) = 0$$

Solving the equation,

>>> solve(2*(-2*x**2 + 1)/(2*x + 1) + (-x**2 + 1)/(x + 1))
[1/12 + sqrt(73)/12, -sqrt(73)/12 + 1/12]
$\endgroup$
4
  • 1
    $\begingroup$ +1: since the JNF is not computable, it suffices to use "Schur decomposition" which is much easier to compute --- the rest of the argument applies unchanged. (link: en.wikipedia.org/wiki/Schur_decomposition ) $\endgroup$ – Suvrit Jan 18 '17 at 19:25
  • 1
    $\begingroup$ @Suvrit Why "not computable"? What does that mean? I have computed the Jordan decomposition using paper & pencil. When I learned it, I was told that it was extremely sensitive to perturbations. Examples convinced me of that claim. $\endgroup$ – Rodrigo de Azevedo Jan 18 '17 at 20:21
  • 1
    $\begingroup$ Yes, what I meant by not computable is that JNF is notorious for being sensitive to numerical concerns, and hence I would not embed it as a subroutine anywhere; as noted in my comment, using your observation, we might as well use Schur decomposition, which is easy to compute and also works, so not only do we not care about diagonalizable, we also do not care about JNF here :-) $\endgroup$ – Suvrit Jan 18 '17 at 22:42
  • 1
    $\begingroup$ @Suvrit I edited my answer and used the Schur decomposition instead, as you recommended. $\endgroup$ – Rodrigo de Azevedo Jan 19 '17 at 10:22

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.