Calculating log2 of a matrix

Question

I'm trying to calculate the base-2-logarithm (in the sense of the inverse function of the matrix-exponential with base 2 as described here, not the elementwise one) of a matrix with python. As log() takes the element-wise logarithm, I did some Google research to find something that is applicable for my task. My research only gave me the scipy function logm(X), which gives back the natural matrix logarithm of the matrix X.

As mentioned above I need to find the base-2-logarithm of a matrix with python. Of course I know the formula $log_a(x) = ln(x)/ln(a)$ where ln is the natural logarithm, but as far as I understand it, this only works for scalar arguments X (Correct me if I'm wrong). At least I haven't seen any argument why this should hold for matrices too.

So does anybody know whether there exists such a built in matrix-log2 function?

Alternatively: Since I've worked a little with Mathematica a few years ago, I know the function MatrixFunction[], which might be a step towards a solution of my problem (as discussed here) and now I'm wondering whether such a function does exist in Python too?

Thanks for your help!

PS: In my opinion the solution presented here doesn't work

Answer 1

Why do you think it works only for scalar argument? Let's play with simple rotation matrix

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

def log2M(a): # base 2 matrix logarithm
    return logm(a)/np.log(2.0)

def p2M(a): # base 2 matrix exponent
    return expm(np.log(2.0)*a)


alpha = 0.6
a = np.array([[np.cos(alpha), -np.sin(alpha)], [np.sin(alpha), np.cos(alpha)]])
print(a)
q = log2M(a)
print(q)

Here you could see it produced reasonable output which looks like

[      0      -alpha/log(2)]
[alpha/log(2)    0         ]

you could compare it with wikipedia

And we could compute it back as 1. Scaled by log(2) exponent 2. Manual by Taylor expansion

And both methods prints the same output, code below

f = 1.0
r = np.array([[1.0, 0.0], [0.0, 1.0]])
eq = np.array([[1.0, 0.0], [0.0, 1.0]])
for k in range(1, 10):
    r = np.dot(r, q)
    f = f * np.float64(k) / np.log(2.0)
    eq = eq + r / f

print(eq)
print(p2M(q))

Output

[[ 0.82533562 -0.56464247]
 [ 0.56464247  0.82533562]]
[[ 0.82533561 -0.56464247]
 [ 0.56464247  0.82533561]]

Everything looks consistent to me