How can I efficiently find an anti-symmetric generator of a special orthogonal matrix?

Question

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.

Answer 1

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$.