diff options
author | Matti Picus <matti.picus@gmail.com> | 2020-02-06 11:06:29 +0200 |
---|---|---|
committer | GitHub <noreply@github.com> | 2020-02-06 11:06:29 +0200 |
commit | 4e5882eb8b11caaf90633bd775b930a314f80deb (patch) | |
tree | 9c9144ece513b7ee1216e425cab8dc20ab457f39 | |
parent | 96fa7141022edf10da6a5a34101d1e8f8920b442 (diff) | |
parent | b202aad8d8e138d4c4cb8ccc590b87e1173f45bc (diff) | |
download | numpy-4e5882eb8b11caaf90633bd775b930a314f80deb.tar.gz |
Merge pull request #15468 from hameerabbasi/fix-svd-sorted
BUG: Fix for SVD not always sorted with hermitian=True
-rw-r--r-- | numpy/linalg/linalg.py | 21 | ||||
-rw-r--r-- | numpy/linalg/tests/test_linalg.py | 8 |
2 files changed, 21 insertions, 8 deletions
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 0964e807b..1d572f876 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -24,7 +24,7 @@ from numpy.core import ( add, multiply, sqrt, fastCopyAndTranspose, sum, isfinite, finfo, errstate, geterrobj, moveaxis, amin, amax, product, abs, atleast_2d, intp, asanyarray, object_, matmul, - swapaxes, divide, count_nonzero, isnan, sign + swapaxes, divide, count_nonzero, isnan, sign, argsort, sort ) from numpy.core.multiarray import normalize_axis_index from numpy.core.overrides import set_module @@ -1608,24 +1608,29 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False): True """ + import numpy as _nx a, wrap = _makearray(a) if hermitian: - # note: lapack returns eigenvalues in reverse order to our contract. - # reversing is cheap by design in numpy, so we do so to be consistent + # note: lapack svd returns eigenvalues with s ** 2 sorted descending, + # but eig returns s sorted ascending, so we re-order the eigenvalues + # and related arrays to have the correct order if compute_uv: s, u = eigh(a) - s = s[..., ::-1] - u = u[..., ::-1] - # singular values are unsigned, move the sign into v - vt = transpose(u * sign(s)[..., None, :]).conjugate() + sgn = sign(s) s = abs(s) + sidx = argsort(s)[..., ::-1] + sgn = _nx.take_along_axis(sgn, sidx, axis=-1) + s = _nx.take_along_axis(s, sidx, axis=-1) + u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1) + # singular values are unsigned, move the sign into v + vt = transpose(u * sgn[..., None, :]).conjugate() return wrap(u), s, wrap(vt) else: s = eigvalsh(a) s = s[..., ::-1] s = abs(s) - return s + return sort(s)[..., ::-1] _assert_stacked_2d(a) t, result_t = _commonType(a) diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py index ae72c4a38..04f5c3d19 100644 --- a/numpy/linalg/tests/test_linalg.py +++ b/numpy/linalg/tests/test_linalg.py @@ -680,6 +680,14 @@ class SVDHermitianCases(HermitianTestCase, HermitianGeneralizedTestCase): assert_allclose(a, dot_generalized(np.asarray(u) * np.asarray(s)[..., None, :], np.asarray(vt)), rtol=get_rtol(u.dtype)) + def hermitian(mat): + axes = list(range(mat.ndim)) + axes[-1], axes[-2] = axes[-2], axes[-1] + return np.conj(np.transpose(mat, axes=axes)) + + assert_almost_equal(np.matmul(u, hermitian(u)), np.broadcast_to(np.eye(u.shape[-1]), u.shape)) + assert_almost_equal(np.matmul(vt, hermitian(vt)), np.broadcast_to(np.eye(vt.shape[-1]), vt.shape)) + assert_equal(np.sort(s)[..., ::-1], s) assert_(consistent_subclass(u, a)) assert_(consistent_subclass(vt, a)) |