[SciPy-Dev] scipy.linalg.svd for n-dimensional arrays in stacked format

[Ilhan Polat]

Hi Melissa,

Yes this is a much-sought after feature not only just SVD but for all
solve-alikes and eig-alikes. I was actually waiting for that magical moment
that I would have some time to work on it but kept on procrastinating. So
by all means, go for it, if you feel like it.

Best
ilhan



On Tue, Feb 18, 2020 at 4:39 PM Melissa Mendonça <melissawm at gmail.com>
wrote:

> Hello all,
>
> I just came across this today. While we generally encourage users to use
> scipy.linalg instead of numpy.linalg, I realized that numpy.linalg.svd
> works for general n-dimensional arrays in "stacked format" (see [1, 2])
> while scipy.linalg.svd does not.
>
> Example:
> =========================
> >>> import numpy as np
> >>> from scipy.linalg import svd
> >>> A = np.random.rand(2,3,4)
> >>> U, s, Vt = np.linalg.svd(A) # works fine
> >>> W, sigma, Yt = svd(A)
>
>
> ---------------------------------------------------------------------------
> ValueError                                Traceback (most recent call last)
> <ipython-input-6-fc853a2e096f> in <module>
> ----> 1 W, sigma, Yt = svd(A)
>
> /opt/miniconda/envs/normal/lib/python3.7/site-packages/scipy/linalg/decomp_svd.py
> in svd(a, full_matrices, compute_uv, overwrite_a, check_finite,
> lapack_driver)
>     109     a1 = _asarray_validated(a, check_finite=check_finite)
>     110     if len(a1.shape) != 2:
> --> 111         raise ValueError('expected matrix')
>     112     m, n = a1.shape
>     113     overwrite_a = overwrite_a or (_datacopied(a1, a))
>
> ValueError: expected matrix
> =========================
>
> My question is if this is a design decision or a missing feature. If this
> is an interesting funcionality that should be added to SciPy, I can look
> into it.
>
> Thanks,
>
> --
> Melissa Weber Mendonça
>
> [1] In the docs for numpy.linalg.svd:  "If `a` has more than two
> dimensions, then broadcasting rules apply, as explained in
> :ref:`routines.linalg-broadcasting`. This means that SVD is working in
> "stacked" mode: it iterates over all indices of the first ``a.ndim - 2``
> dimensions and for each combination SVD is applied to the last two indices.
> The matrix `a` can be reconstructed from the decomposition with either ``(u
> * s[..., None, :]) @ vh`` or ``u @ (s[..., None] * vh)``. (The ``@``
> operator can be replaced by the function ``np.matmul`` for python versions
> below 3.5.)
> [2] https://numpy.org/devdocs/reference/generated/numpy.linalg.svd.html
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at python.org
> https://mail.python.org/mailman/listinfo/scipy-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20200218/959e7c8e/attachment.html>