Is there a matrix logarithm in Sage?

Question

I just started using SageMath to do some linear algebra computations. The exponential map $\exp$ is built in for matrices. Is the logarithm

$$\log(u) = (u-I) - \frac{1}{2}(u-I)^2 + \frac{1}{3}(u-I)^3 - \cdots$$

also built in? For nilpotent (unipotent) matrices, $\exp$ ($\log$) gives finite sums, and I want to make use of the identities

$$\exp \circ \log(u) = u, \log \circ \exp(X) = X$$

for unipotent $u$ and nilpotent $X$.

Answer 1

You can convert your matrix to numpy and use logm of SciPy:

from scipy.linalg import logm
import numpy as np
M = random_matrix(RDF,4,4)  # your matrix
M_np = np.matrix(M,dtype='float64')
logm(M_np)

To convert numpy array to Sage matrix use:

matrix(np.asmatrix(logm(M_np)))

Obs: convert numpy array to sage matrix directly gave error in some cases, so one should convert to numpy matrix before converting to sage matrix.

If you need log often it may be convenient to do a function:

def log_matrix(M):
    from scipy.linalg import logm
    import numpy as np
    M_np = np.matrix(M,dtype='float64')
    return matrix(np.asmatrix(logm(M_np)))

with Input and Output sage matrix.

Note that exp(log_matrix(M))-M and log_matrix(exp(M))-M differs from the null matrix with an error of the order of $1e-15$.

Answer 2

Following @Travis's suggestion, a one-liner to compute the logarithm of a unipotent matrix.

Define a unipotent matrix:

sage: u = matrix([[1, 1], [0, 1]])
sage: u
[1 1]
[0 1]

Check that it is unipotent:

sage: (u-1)^u.ncols()
[0 0]
[0 0]

Compute its logarithm:

sage: log_u = sum((-1)^(k+1)/k*(u-1)^k for k in range(1, u.ncols()))
sage: log_u
[0 1]
[0 0]

Turn this into a function (with a check that the argument is a unipotent matrix):

def log_of_unipotent(u):
    if (u-1)^u.ncols() != 0:
        raise ValueError('This matrix is not unipotent')
    return sum((-1)^(k+1)/k*(u-1)^k for k in range(1, u.ncols()))

Try it out taking inspiration from the Wikipedia page on nilpotent matrices:

sage: n = matrix([[5, -3, 2], [15, -9, 6], [10, -6, 4]])
sage: n
[ 5 -3  2]
[15 -9  6]
[10 -6  4]
sage: n^n.ncols() == 0  # check n is nilpotent
True
sage: u = n.exp()
sage: u
[ 6 -3  2]
[15 -8  6]
[10 -6  5]
sage: (u-1)^u.ncols() == 0  # check u is unipotent
True
sage: log_of_unipotent(exp(n)) == n
True
sage: exp(log_of_unipotent(u)) == u
True

Note that we stick here to unipotent matrices, as suggested by the question. For more general considerations see for example the Wikipedia page on the logarithm of a matrix.

Finally, depending on one's taste, one could

• replace (-1)^(k+1) by (1 if k%2 else -1)
• replace (u - 1) by (u - u^0) or (u - u.parent().identity_matrix())