diff options
Diffstat (limited to 'numpy/linalg/linalg.py')
| -rw-r--r-- | numpy/linalg/linalg.py | 203 |
1 files changed, 116 insertions, 87 deletions
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 8363d7377..17e84be2d 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -26,7 +26,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 + swapaxes, divide, count_nonzero, isnan, sign ) from numpy.core.multiarray import normalize_axis_index from numpy.core.overrides import set_module @@ -377,7 +377,7 @@ def solve(a, b): >>> b = np.array([9,8]) >>> x = np.linalg.solve(a, b) >>> x - array([ 2., 3.]) + array([2., 3.]) Check that the solution is correct: @@ -535,10 +535,10 @@ def inv(a): >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]]) >>> inv(a) - array([[[-2. , 1. ], - [ 1.5, -0.5]], - [[-5. , 2. ], - [ 3. , -1. ]]]) + array([[[-2. , 1. ], + [ 1.5 , -0.5 ]], + [[-1.25, 0.75], + [ 0.75, -0.25]]]) """ a, wrap = _makearray(a) @@ -730,21 +730,21 @@ def cholesky(a): -------- >>> A = np.array([[1,-2j],[2j,5]]) >>> A - array([[ 1.+0.j, 0.-2.j], + array([[ 1.+0.j, -0.-2.j], [ 0.+2.j, 5.+0.j]]) >>> L = np.linalg.cholesky(A) >>> L - array([[ 1.+0.j, 0.+0.j], - [ 0.+2.j, 1.+0.j]]) + array([[1.+0.j, 0.+0.j], + [0.+2.j, 1.+0.j]]) >>> np.dot(L, L.T.conj()) # verify that L * L.H = A - array([[ 1.+0.j, 0.-2.j], - [ 0.+2.j, 5.+0.j]]) + array([[1.+0.j, 0.-2.j], + [0.+2.j, 5.+0.j]]) >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like? >>> np.linalg.cholesky(A) # an ndarray object is returned - array([[ 1.+0.j, 0.+0.j], - [ 0.+2.j, 1.+0.j]]) + array([[1.+0.j, 0.+0.j], + [0.+2.j, 1.+0.j]]) >>> # But a matrix object is returned if A is a matrix object - >>> LA.cholesky(np.matrix(A)) + >>> np.linalg.cholesky(np.matrix(A)) matrix([[ 1.+0.j, 0.+0.j], [ 0.+2.j, 1.+0.j]]) @@ -878,9 +878,9 @@ def qr(a, mode='reduced'): [1, 1], [2, 1]]) >>> b = np.array([1, 0, 2, 1]) - >>> q, r = LA.qr(A) + >>> q, r = np.linalg.qr(A) >>> p = np.dot(q.T, b) - >>> np.dot(LA.inv(r), p) + >>> np.dot(np.linalg.inv(r), p) array([ 1.1e-16, 1.0e+00]) """ @@ -1049,7 +1049,7 @@ def eigvals(a): >>> A = np.dot(Q, D) >>> A = np.dot(A, Q.T) >>> LA.eigvals(A) - array([ 1., -1.]) + array([ 1., -1.]) # random """ a, wrap = _makearray(a) @@ -1131,24 +1131,24 @@ def eigvalsh(a, UPLO='L'): >>> from numpy import linalg as LA >>> a = np.array([[1, -2j], [2j, 5]]) >>> LA.eigvalsh(a) - array([ 0.17157288, 5.82842712]) + array([ 0.17157288, 5.82842712]) # may vary >>> # demonstrate the treatment of the imaginary part of the diagonal >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]]) >>> a - array([[ 5.+2.j, 9.-2.j], - [ 0.+2.j, 2.-1.j]]) + array([[5.+2.j, 9.-2.j], + [0.+2.j, 2.-1.j]]) >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals() >>> # with: >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]]) >>> b - array([[ 5.+0.j, 0.-2.j], - [ 0.+2.j, 2.+0.j]]) + array([[5.+0.j, 0.-2.j], + [0.+2.j, 2.+0.j]]) >>> wa = LA.eigvalsh(a) >>> wb = LA.eigvals(b) >>> wa; wb - array([ 1., 6.]) - array([ 6.+0.j, 1.+0.j]) + array([1., 6.]) + array([6.+0.j, 1.+0.j]) """ UPLO = UPLO.upper() @@ -1264,19 +1264,19 @@ def eig(a): >>> w, v = LA.eig(np.diag((1, 2, 3))) >>> w; v - array([ 1., 2., 3.]) - array([[ 1., 0., 0.], - [ 0., 1., 0.], - [ 0., 0., 1.]]) + array([1., 2., 3.]) + array([[1., 0., 0.], + [0., 1., 0.], + [0., 0., 1.]]) Real matrix possessing complex e-values and e-vectors; note that the e-values are complex conjugates of each other. >>> w, v = LA.eig(np.array([[1, -1], [1, 1]])) >>> w; v - array([ 1. + 1.j, 1. - 1.j]) - array([[ 0.70710678+0.j , 0.70710678+0.j ], - [ 0.00000000-0.70710678j, 0.00000000+0.70710678j]]) + array([1.+1.j, 1.-1.j]) + array([[0.70710678+0.j , 0.70710678-0.j ], + [0. -0.70710678j, 0. +0.70710678j]]) Complex-valued matrix with real e-values (but complex-valued e-vectors); note that a.conj().T = a, i.e., a is Hermitian. @@ -1284,9 +1284,9 @@ def eig(a): >>> a = np.array([[1, 1j], [-1j, 1]]) >>> w, v = LA.eig(a) >>> w; v - array([ 2.00000000e+00+0.j, 5.98651912e-36+0.j]) # i.e., {2, 0} - array([[ 0.00000000+0.70710678j, 0.70710678+0.j ], - [ 0.70710678+0.j , 0.00000000+0.70710678j]]) + array([2.+0.j, 0.+0.j]) + array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary + [ 0.70710678+0.j , -0. +0.70710678j]]) Be careful about round-off error! @@ -1294,9 +1294,9 @@ def eig(a): >>> # Theor. e-values are 1 +/- 1e-9 >>> w, v = LA.eig(a) >>> w; v - array([ 1., 1.]) - array([[ 1., 0.], - [ 0., 1.]]) + array([1., 1.]) + array([[1., 0.], + [0., 1.]]) """ a, wrap = _makearray(a) @@ -1392,49 +1392,49 @@ def eigh(a, UPLO='L'): >>> from numpy import linalg as LA >>> a = np.array([[1, -2j], [2j, 5]]) >>> a - array([[ 1.+0.j, 0.-2.j], + array([[ 1.+0.j, -0.-2.j], [ 0.+2.j, 5.+0.j]]) >>> w, v = LA.eigh(a) >>> w; v - array([ 0.17157288, 5.82842712]) - array([[-0.92387953+0.j , -0.38268343+0.j ], - [ 0.00000000+0.38268343j, 0.00000000-0.92387953j]]) + array([0.17157288, 5.82842712]) + array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary + [ 0. +0.38268343j, 0. -0.92387953j]]) >>> np.dot(a, v[:, 0]) - w[0] * v[:, 0] # verify 1st e-val/vec pair - array([2.77555756e-17 + 0.j, 0. + 1.38777878e-16j]) + array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j]) >>> np.dot(a, v[:, 1]) - w[1] * v[:, 1] # verify 2nd e-val/vec pair - array([ 0.+0.j, 0.+0.j]) + array([0.+0.j, 0.+0.j]) >>> A = np.matrix(a) # what happens if input is a matrix object >>> A - matrix([[ 1.+0.j, 0.-2.j], + matrix([[ 1.+0.j, -0.-2.j], [ 0.+2.j, 5.+0.j]]) >>> w, v = LA.eigh(A) >>> w; v - array([ 0.17157288, 5.82842712]) - matrix([[-0.92387953+0.j , -0.38268343+0.j ], - [ 0.00000000+0.38268343j, 0.00000000-0.92387953j]]) + array([0.17157288, 5.82842712]) + matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary + [ 0. +0.38268343j, 0. -0.92387953j]]) >>> # demonstrate the treatment of the imaginary part of the diagonal >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]]) >>> a - array([[ 5.+2.j, 9.-2.j], - [ 0.+2.j, 2.-1.j]]) + array([[5.+2.j, 9.-2.j], + [0.+2.j, 2.-1.j]]) >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with: >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]]) >>> b - array([[ 5.+0.j, 0.-2.j], - [ 0.+2.j, 2.+0.j]]) + array([[5.+0.j, 0.-2.j], + [0.+2.j, 2.+0.j]]) >>> wa, va = LA.eigh(a) >>> wb, vb = LA.eig(b) >>> wa; wb - array([ 1., 6.]) - array([ 6.+0.j, 1.+0.j]) + array([1., 6.]) + array([6.+0.j, 1.+0.j]) >>> va; vb - array([[-0.44721360-0.j , -0.89442719+0.j ], - [ 0.00000000+0.89442719j, 0.00000000-0.4472136j ]]) - array([[ 0.89442719+0.j , 0.00000000-0.4472136j], - [ 0.00000000-0.4472136j, 0.89442719+0.j ]]) + array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary + [ 0. +0.89442719j, 0. -0.4472136j ]]) + array([[ 0.89442719+0.j , -0. +0.4472136j], + [-0. +0.4472136j, 0.89442719+0.j ]]) """ UPLO = UPLO.upper() if UPLO not in ('L', 'U'): @@ -1461,12 +1461,12 @@ def eigh(a, UPLO='L'): # Singular value decomposition -def _svd_dispatcher(a, full_matrices=None, compute_uv=None): +def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None): return (a,) @array_function_dispatch(_svd_dispatcher) -def svd(a, full_matrices=True, compute_uv=True): +def svd(a, full_matrices=True, compute_uv=True, hermitian=False): """ Singular Value Decomposition. @@ -1504,6 +1504,12 @@ def svd(a, full_matrices=True, compute_uv=True): size as those of the input `a`. The size of the last two dimensions depends on the value of `full_matrices`. Only returned when `compute_uv` is True. + hermitian : bool, optional + If True, `a` is assumed to be Hermitian (symmetric if real-valued), + enabling a more efficient method for finding singular values. + Defaults to False. + + ..versionadded:: 1.17.0 Raises ------ @@ -1590,6 +1596,24 @@ def svd(a, full_matrices=True, compute_uv=True): """ 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 + 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() + s = abs(s) + return wrap(u), s, wrap(vt) + else: + s = eigvalsh(a) + s = s[..., ::-1] + s = abs(s) + return s + _assertRankAtLeast2(a) t, result_t = _commonType(a) @@ -1705,9 +1729,9 @@ def cond(x, p=None): >>> LA.cond(a, 2) 1.4142135623730951 >>> LA.cond(a, -2) - 0.70710678118654746 + 0.70710678118654746 # may vary >>> min(LA.svd(a, compute_uv=0))*min(LA.svd(LA.inv(a), compute_uv=0)) - 0.70710678118654746 + 0.70710678118654746 # may vary """ x = asarray(x) # in case we have a matrix @@ -1844,10 +1868,7 @@ def matrix_rank(M, tol=None, hermitian=False): M = asarray(M) if M.ndim < 2: return int(not all(M==0)) - if hermitian: - S = abs(eigvalsh(M)) - else: - S = svd(M, compute_uv=False) + S = svd(M, compute_uv=False, hermitian=hermitian) if tol is None: tol = S.max(axis=-1, keepdims=True) * max(M.shape[-2:]) * finfo(S.dtype).eps else: @@ -1857,12 +1878,12 @@ def matrix_rank(M, tol=None, hermitian=False): # Generalized inverse -def _pinv_dispatcher(a, rcond=None): +def _pinv_dispatcher(a, rcond=None, hermitian=None): return (a,) @array_function_dispatch(_pinv_dispatcher) -def pinv(a, rcond=1e-15): +def pinv(a, rcond=1e-15, hermitian=False): """ Compute the (Moore-Penrose) pseudo-inverse of a matrix. @@ -1882,6 +1903,12 @@ def pinv(a, rcond=1e-15): Singular values smaller (in modulus) than `rcond` * largest_singular_value (again, in modulus) are set to zero. Broadcasts against the stack of matrices + hermitian : bool, optional + If True, `a` is assumed to be Hermitian (symmetric if real-valued), + enabling a more efficient method for finding singular values. + Defaults to False. + + ..versionadded:: 1.17.0 Returns ------- @@ -1935,7 +1962,7 @@ def pinv(a, rcond=1e-15): res = empty(a.shape[:-2] + (n, m), dtype=a.dtype) return wrap(res) a = a.conjugate() - u, s, vt = svd(a, full_matrices=False) + u, s, vt = svd(a, full_matrices=False, hermitian=hermitian) # discard small singular values cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True) @@ -2002,7 +2029,7 @@ def slogdet(a): >>> a = np.array([[1, 2], [3, 4]]) >>> (sign, logdet) = np.linalg.slogdet(a) >>> (sign, logdet) - (-1, 0.69314718055994529) + (-1, 0.69314718055994529) # may vary >>> sign * np.exp(logdet) -2.0 @@ -2074,7 +2101,7 @@ def det(a): >>> a = np.array([[1, 2], [3, 4]]) >>> np.linalg.det(a) - -2.0 + -2.0 # may vary Computing determinants for a stack of matrices: @@ -2181,15 +2208,15 @@ def lstsq(a, b, rcond="warn"): [ 3., 1.]]) >>> m, c = np.linalg.lstsq(A, y, rcond=None)[0] - >>> print(m, c) - 1.0 -0.95 + >>> m, c + (1.0 -0.95) # may vary Plot the data along with the fitted line: >>> import matplotlib.pyplot as plt - >>> plt.plot(x, y, 'o', label='Original data', markersize=10) - >>> plt.plot(x, m*x + c, 'r', label='Fitted line') - >>> plt.legend() + >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10) + >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line') + >>> _ = plt.legend() >>> plt.show() """ @@ -2367,7 +2394,7 @@ def norm(x, ord=None, axis=None, keepdims=False): >>> from numpy import linalg as LA >>> a = np.arange(9) - 4 >>> a - array([-4, -3, -2, -1, 0, 1, 2, 3, 4]) + array([-4, -3, -2, ..., 2, 3, 4]) >>> b = a.reshape((3, 3)) >>> b array([[-4, -3, -2], @@ -2403,13 +2430,13 @@ def norm(x, ord=None, axis=None, keepdims=False): 7.3484692283495345 >>> LA.norm(a, -2) - nan + 0.0 >>> LA.norm(b, -2) - 1.8570331885190563e-016 + 1.8570331885190563e-016 # may vary >>> LA.norm(a, 3) - 5.8480354764257312 + 5.8480354764257312 # may vary >>> LA.norm(a, -3) - nan + 0.0 Using the `axis` argument to compute vector norms: @@ -2584,18 +2611,20 @@ def multi_dot(arrays): >>> from numpy.linalg import multi_dot >>> # Prepare some data - >>> A = np.random.random(10000, 100) - >>> B = np.random.random(100, 1000) - >>> C = np.random.random(1000, 5) - >>> D = np.random.random(5, 333) + >>> A = np.random.random((10000, 100)) + >>> B = np.random.random((100, 1000)) + >>> C = np.random.random((1000, 5)) + >>> D = np.random.random((5, 333)) >>> # the actual dot multiplication - >>> multi_dot([A, B, C, D]) + >>> _ = multi_dot([A, B, C, D]) instead of:: - >>> np.dot(np.dot(np.dot(A, B), C), D) + >>> _ = np.dot(np.dot(np.dot(A, B), C), D) + ... >>> # or - >>> A.dot(B).dot(C).dot(D) + >>> _ = A.dot(B).dot(C).dot(D) + ... Notes ----- |