exp not working on SO(3)

[mateuszbaran]

I was trying to compute the exponential map on SO(3) but I couldn't get it to work in any example (using the numpy backend). Here is my code:

import pymanopt, numpy

manifold = pymanopt.manifolds.special_orthogonal_group.SpecialOrthogonalGroup(3)
p = manifold.random_point()
X = manifold.random_tangent_vector(p)
manifold.exp(p, X)

and that results in the following error message:

Traceback (most recent call last):
  File "<input>", line 5, in <module>
  File "/home/mateusz/repos/benchmark_tests/venv/lib/python3.8/site-packages/pymanopt/manifolds/special_orthogonal_group.py", line 123, in exp
    return point @ multiexpm(tangent_vector)
  File "/home/mateusz/repos/benchmark_tests/venv/lib/python3.8/site-packages/pymanopt/tools/multi.py", line 65, in multiexpm
    return np.vectorize(scipy.linalg.expm, signature="(m,m)->(m,m)")(A)
  File "/home/mateusz/repos/benchmark_tests/venv/lib/python3.8/site-packages/numpy/lib/function_base.py", line 2328, in __call__
    return self._vectorize_call(func=func, args=vargs)
  File "/home/mateusz/repos/benchmark_tests/venv/lib/python3.8/site-packages/numpy/lib/function_base.py", line 2402, in _vectorize_call
    res = self._vectorize_call_with_signature(func, args)
  File "/home/mateusz/repos/benchmark_tests/venv/lib/python3.8/site-packages/numpy/lib/function_base.py", line 2442, in _vectorize_call_with_signature
    results = func(*(arg[index] for arg in args))
  File "/home/mateusz/repos/benchmark_tests/venv/lib/python3.8/site-packages/scipy/linalg/_matfuncs.py", line 335, in expm
    lu = bandwidth(aw)
  File "_cythonized_array_utils.pyx", line 104, in scipy.linalg._cythonized_array_utils.bandwidth
  File "_cythonized_array_utils.pyx", line 114, in scipy.linalg._cythonized_array_utils.bandwidth_c
  File "stringsource", line 660, in View.MemoryView.memoryview_cwrapper
  File "stringsource", line 350, in View.MemoryView.memoryview.__cinit__
ValueError: buffer source array is read-only

How could I fix that?

[qrebjock]

Hello mateuszbaran,

I don't have any error when I run your code.

What versions of numpy, scipy and cython do you have?
Mine are 1.19.4, 1.5.4 and 0.29.21 respectively.

Did you also pull the last version of pymanopt?

[mateuszbaran]

I'm using CPython 3.8, numpy 1.19.5, scipy 1.9.0 and no cython as far as I can tell. I'm using pymanopt 2.0.0 installed in a virtual environment using pip. I have downgraded scipy to 1.5.4 and it works, so it looks like scipy is causing this.

[nkoep]

I can reproduce this as well on with scipy 1.9.0. Downgrading to 1.8.1 solves the issue. Apparently scipy 1.9.0 added support for calling expm on stacked matrices (see scipy/scipy#15079 and https://scipy.github.io/devdocs/release.1.9.0.html#scipy-linalg-improvements). So for an array of shape (...,n,n) with k >= 2 dimensions, they apply expm to each of the k-2 dimenions automatically. In older versions we faked this behavior by wrapping scipy.linalg.expm with np.vectorize. I'll raise an issue in the scipy repo to check whether this exception is expected when vectorizing such a function.

[nkoep]

Here's the ticket in case anyone wants to track it: scipy/scipy#16792.

[nkoep]

I just released pymanopt 2.0.1 on pypi which should fix this issue. Let's see if the issue gets fixed in scipy upstream so we can remove our workaround again, but for now this should do. Thanks for raising this!

[mateuszbaran]

Thank you for looking into the issue!