4
$\begingroup$

Given a special orthogonal matrix $O$ (i.e: $OO^T = 1$ and $\det(O) = 1$), I am trying to efficiently find a matrix $X$ such that $O = e^X$ and $X = -X^T$ using Python (NumPy & SciPy).

One obvious way to do this is to calculate the logarithm of $O$ using the function scipy.linalg.logm, and indeed this seems to work for almost all matrices. Namely, I verified that this works by running the following code

from scipy.stats import ortho_group
from scipy.linalg import expm, logm
import numpy as np

dim = 16

while True:
    # generate random SO(dim) matrix
    O = ortho_group.rvs(dim)
    if not np.abs(np.linalg.det(O) - 1) < 1e-8:
        O[:,0] = -O[:,0]
    
    # compute generator and check if it is anti-symmetric
    X = logm(O)

    assert np.allclose(X, -X.T)

However, one specific matrix I am trying to find the generator of is the following:

O = np.array([[-4.89378834e-01, -1.11359943e-15, -8.72071303e-01,
         2.08652704e-15],
       [-2.22044605e-16,  5.40669135e-01, -2.55351296e-15,
        -8.41235334e-01],
       [-8.72071303e-01, -3.33066907e-16,  4.89378834e-01,
        -1.55431223e-15],
       [ 8.79435558e-16, -8.41235334e-01, -7.12898969e-16,
        -5.40669135e-01]])

and if I compute X = logm(O), I get that np.allclose(-X.T, X) returns False, whereas np.allclose(-X.conj().T, X) returns True. Thus, I am finding an anti-Hermitian generator, rather than an anti-symmetric generator. Unfortunately this is not sufficient for my overall application. This motivates the question: is there a simple, efficient way to find an anti-symmetric generator of a special orthogonal matrix?

I should point out here: of course, I could in principle do this via non-linear optimization of the anti-symmetric generator using scipy.optimize.minimize on the cost function $f(X) = ||O - e^X||$. However, this is quite slow, and it's not guaranteed that I will get a numerically exact solution.

$\endgroup$
4
  • $\begingroup$ I modified your infinite loop in your code and ran it several times (10000 iterations in each case) and all the times it gave matrices that satisfied your criteria. One difference between your example and the other matrices is the eigenvalues; your 4 by 4 matrix has purely real eigenvalues. $\endgroup$
    – nicoguaro
    Commented Jul 2 at 17:26
  • $\begingroup$ I think your problem is negative eigenvalues. What do you want your function to return on the matrix [-1]? $\endgroup$ Commented Jul 2 at 21:08
  • $\begingroup$ @FedericoPoloni I would not call my function on the matrix [-1], since I am checking beforehand that $det(O) = 1$ for any matrix I call the function on. If this weren't the case I would scale a column by $-1$ to make it so. Mathematically I wouldn't be able to find an anti-symmetric generator for O = [-1]. $\endgroup$ Commented Jul 3 at 0:23
  • $\begingroup$ @nicoguaro That's a good point. It might be a problem that only arises when $O$ has real eigenvalues in $\{-1,1\}$. That would explain why the loop check succeeds with probability $\approx 1$. $\endgroup$ Commented Jul 3 at 0:24

1 Answer 1

7
$\begingroup$

As I had noted in comments (now deleted), the problem is when you have eigenvalues equal to $-1$; those must come in pairs, since $\det O = 1$. If $O = \begin{bmatrix}-1 \\ & -1\end{bmatrix}$, then you want to return $X = \begin{bmatrix}0 & \pi \\ -\pi & 0\end{bmatrix}$, while logm returns $X = \begin{bmatrix}-\pi i \\ & -\pi i\end{bmatrix}$, a different branch of the matrix logarithm.

Returning the first branch is trickier, because it's not a "primary function" in matrix function terms, i.e., $X$ cannot be written as a polynomial in $O$. Any method to compute logm based on a polynomial or rational approximation of the logarithm won't work. Luckily, though, $O$ is normal, so it can be diagonalized with a perfectly stable orthogonal change of basis. In practice, you can work as follows.

  1. Start by computing a real Schur form $O = UTU^\top$. Since $O$ is orthogonal, and hence normal, $T$ is going to be block diagonal, with orthogonal $1\times 1$ or $2\times 2$ blocks in $T$.

  2. We want to compute the correct branch of $\log T$, block by block, so that it is real antisymmetric. On the $2\times 2$ blocks, the function logm already returns what you need. On the $1\times 1$ block [1], just return $0$. For the $1\times 1$ blocks [-1], you need to combine them in pairs and return the branch $$ \log \begin{bmatrix}-1\\ & -1\end{bmatrix} = \begin{bmatrix}&\pi\\ -\pi\end{bmatrix} $$ rather than what logm would return. To form these pairs, reorder the columns of $U$ as well as the rows and columns of $T$, so that the $-1$ are all adjacent.

  3. Once you have $\log T$, you can compute the product $\log O = U (\log T) U^\top$.

$\endgroup$
3
  • $\begingroup$ The MATLAB command schur() will find the Schur decomposition of your matrix. $\endgroup$ Commented Jul 5 at 17:35
  • 3
    $\begingroup$ In scipy, there's scipy.linalg.schur() $\endgroup$ Commented Jul 5 at 18:09
  • 2
    $\begingroup$ That's a great answer! Thanks! $\endgroup$ Commented Jul 5 at 20:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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