1"""Lite version of scipy.linalg.
2
3Notes
4-----
5This module is a lite version of the linalg.py module in SciPy which
6contains high-level Python interface to the LAPACK library. The lite
7version only accesses the following LAPACK functions: dgesv, zgesv,
8dgeev, zgeev, dgesdd, zgesdd, dgelsd, zgelsd, dsyevd, zheevd, dgetrf,
9zgetrf, dpotrf, zpotrf, dgeqrf, zgeqrf, zungqr, dorgqr.
10"""
11
12__all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv',
13 'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'slogdet', 'det',
14 'svd', 'eig', 'eigh', 'lstsq', 'norm', 'qr', 'cond', 'matrix_rank',
15 'LinAlgError', 'multi_dot']
16
17import functools
18import operator
19import warnings
20
21from numpy.core import (
22 array, asarray, zeros, empty, empty_like, intc, single, double,
23 csingle, cdouble, inexact, complexfloating, newaxis, all, Inf, dot,
24 add, multiply, sqrt, sum, isfinite,
25 finfo, errstate, geterrobj, moveaxis, amin, amax, product, abs,
26 atleast_2d, intp, asanyarray, object_, matmul,
27 swapaxes, divide, count_nonzero, isnan, sign, argsort, sort,
28 reciprocal
29)
30from numpy.core.multiarray import normalize_axis_index
31from numpy.core.overrides import set_module
32from numpy.core import overrides
33from numpy.lib.twodim_base import triu, eye
34from numpy.linalg import _umath_linalg
35
36
37array_function_dispatch = functools.partial(
38 overrides.array_function_dispatch, module='numpy.linalg')
39
40
41fortran_int = intc
42
43
44@set_module('numpy.linalg')
45class LinAlgError(Exception):
46 """
47 Generic Python-exception-derived object raised by linalg functions.
48
49 General purpose exception class, derived from Python's exception.Exception
50 class, programmatically raised in linalg functions when a Linear
51 Algebra-related condition would prevent further correct execution of the
52 function.
53
54 Parameters
55 ----------
56 None
57
58 Examples
59 --------
60 >>> from numpy import linalg as LA
61 >>> LA.inv(np.zeros((2,2)))
62 Traceback (most recent call last):
63 File "<stdin>", line 1, in <module>
64 File "...linalg.py", line 350,
65 in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
66 File "...linalg.py", line 249,
67 in solve
68 raise LinAlgError('Singular matrix')
69 numpy.linalg.LinAlgError: Singular matrix
70
71 """
72
73
74def _determine_error_states():
75 errobj = geterrobj()
76 bufsize = errobj[0]
77
78 with errstate(invalid='call', over='ignore',
79 divide='ignore', under='ignore'):
80 invalid_call_errmask = geterrobj()[1]
81
82 return [bufsize, invalid_call_errmask, None]
83
84# Dealing with errors in _umath_linalg
85_linalg_error_extobj = _determine_error_states()
86del _determine_error_states
87
88def _raise_linalgerror_singular(err, flag):
89 raise LinAlgError("Singular matrix")
90
91def _raise_linalgerror_nonposdef(err, flag):
92 raise LinAlgError("Matrix is not positive definite")
93
94def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
95 raise LinAlgError("Eigenvalues did not converge")
96
97def _raise_linalgerror_svd_nonconvergence(err, flag):
98 raise LinAlgError("SVD did not converge")
99
100def _raise_linalgerror_lstsq(err, flag):
101 raise LinAlgError("SVD did not converge in Linear Least Squares")
102
103def _raise_linalgerror_qr(err, flag):
104 raise LinAlgError("Incorrect argument found while performing "
105 "QR factorization")
106
107def get_linalg_error_extobj(callback):
108 extobj = list(_linalg_error_extobj) # make a copy
109 extobj[2] = callback
110 return extobj
111
112def _makearray(a):
113 new = asarray(a)
114 wrap = getattr(a, "__array_prepare__", new.__array_wrap__)
115 return new, wrap
116
117def isComplexType(t):
118 return issubclass(t, complexfloating)
119
120_real_types_map = {single : single,
121 double : double,
122 csingle : single,
123 cdouble : double}
124
125_complex_types_map = {single : csingle,
126 double : cdouble,
127 csingle : csingle,
128 cdouble : cdouble}
129
130def _realType(t, default=double):
131 return _real_types_map.get(t, default)
132
133def _complexType(t, default=cdouble):
134 return _complex_types_map.get(t, default)
135
136def _commonType(*arrays):
137 # in lite version, use higher precision (always double or cdouble)
138 result_type = single
139 is_complex = False
140 for a in arrays:
141 if issubclass(a.dtype.type, inexact):
142 if isComplexType(a.dtype.type):
143 is_complex = True
144 rt = _realType(a.dtype.type, default=None)
145 if rt is None:
146 # unsupported inexact scalar
147 raise TypeError("array type %s is unsupported in linalg" %
148 (a.dtype.name,))
149 else:
150 rt = double
151 if rt is double:
152 result_type = double
153 if is_complex:
154 t = cdouble
155 result_type = _complex_types_map[result_type]
156 else:
157 t = double
158 return t, result_type
159
160
161def _to_native_byte_order(*arrays):
162 ret = []
163 for arr in arrays:
164 if arr.dtype.byteorder not in ('=', '|'):
165 ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('=')))
166 else:
167 ret.append(arr)
168 if len(ret) == 1:
169 return ret[0]
170 else:
171 return ret
172
173
174def _assert_2d(*arrays):
175 for a in arrays:
176 if a.ndim != 2:
177 raise LinAlgError('%d-dimensional array given. Array must be '
178 'two-dimensional' % a.ndim)
179
180def _assert_stacked_2d(*arrays):
181 for a in arrays:
182 if a.ndim < 2:
183 raise LinAlgError('%d-dimensional array given. Array must be '
184 'at least two-dimensional' % a.ndim)
185
186def _assert_stacked_square(*arrays):
187 for a in arrays:
188 m, n = a.shape[-2:]
189 if m != n:
190 raise LinAlgError('Last 2 dimensions of the array must be square')
191
192def _assert_finite(*arrays):
193 for a in arrays:
194 if not isfinite(a).all():
195 raise LinAlgError("Array must not contain infs or NaNs")
196
197def _is_empty_2d(arr):
198 # check size first for efficiency
199 return arr.size == 0 and product(arr.shape[-2:]) == 0
200
201
202def transpose(a):
203 """
204 Transpose each matrix in a stack of matrices.
205
206 Unlike np.transpose, this only swaps the last two axes, rather than all of
207 them
208
209 Parameters
210 ----------
211 a : (...,M,N) array_like
212
213 Returns
214 -------
215 aT : (...,N,M) ndarray
216 """
217 return swapaxes(a, -1, -2)
218
219# Linear equations
220
221def _tensorsolve_dispatcher(a, b, axes=None):
222 return (a, b)
223
224
225@array_function_dispatch(_tensorsolve_dispatcher)
226def tensorsolve(a, b, axes=None):
227 """
228 Solve the tensor equation ``a x = b`` for x.
229
230 It is assumed that all indices of `x` are summed over in the product,
231 together with the rightmost indices of `a`, as is done in, for example,
232 ``tensordot(a, x, axes=x.ndim)``.
233
234 Parameters
235 ----------
236 a : array_like
237 Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals
238 the shape of that sub-tensor of `a` consisting of the appropriate
239 number of its rightmost indices, and must be such that
240 ``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be
241 'square').
242 b : array_like
243 Right-hand tensor, which can be of any shape.
244 axes : tuple of ints, optional
245 Axes in `a` to reorder to the right, before inversion.
246 If None (default), no reordering is done.
247
248 Returns
249 -------
250 x : ndarray, shape Q
251
252 Raises
253 ------
254 LinAlgError
255 If `a` is singular or not 'square' (in the above sense).
256
257 See Also
258 --------
259 numpy.tensordot, tensorinv, numpy.einsum
260
261 Examples
262 --------
263 >>> a = np.eye(2*3*4)
264 >>> a.shape = (2*3, 4, 2, 3, 4)
265 >>> b = np.random.randn(2*3, 4)
266 >>> x = np.linalg.tensorsolve(a, b)
267 >>> x.shape
268 (2, 3, 4)
269 >>> np.allclose(np.tensordot(a, x, axes=3), b)
270 True
271
272 """
273 a, wrap = _makearray(a)
274 b = asarray(b)
275 an = a.ndim
276
277 if axes is not None:
278 allaxes = list(range(0, an))
279 for k in axes:
280 allaxes.remove(k)
281 allaxes.insert(an, k)
282 a = a.transpose(allaxes)
283
284 oldshape = a.shape[-(an-b.ndim):]
285 prod = 1
286 for k in oldshape:
287 prod *= k
288
289 if a.size != prod ** 2:
290 raise LinAlgError(
291 "Input arrays must satisfy the requirement \
292 prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])"
293 )
294
295 a = a.reshape(prod, prod)
296 b = b.ravel()
297 res = wrap(solve(a, b))
298 res.shape = oldshape
299 return res
300
301
302def _solve_dispatcher(a, b):
303 return (a, b)
304
305
306@array_function_dispatch(_solve_dispatcher)
307def solve(a, b):
308 """
309 Solve a linear matrix equation, or system of linear scalar equations.
310
311 Computes the "exact" solution, `x`, of the well-determined, i.e., full
312 rank, linear matrix equation `ax = b`.
313
314 Parameters
315 ----------
316 a : (..., M, M) array_like
317 Coefficient matrix.
318 b : {(..., M,), (..., M, K)}, array_like
319 Ordinate or "dependent variable" values.
320
321 Returns
322 -------
323 x : {(..., M,), (..., M, K)} ndarray
324 Solution to the system a x = b. Returned shape is identical to `b`.
325
326 Raises
327 ------
328 LinAlgError
329 If `a` is singular or not square.
330
331 See Also
332 --------
333 scipy.linalg.solve : Similar function in SciPy.
334
335 Notes
336 -----
337
338 .. versionadded:: 1.8.0
339
340 Broadcasting rules apply, see the `numpy.linalg` documentation for
341 details.
342
343 The solutions are computed using LAPACK routine ``_gesv``.
344
345 `a` must be square and of full-rank, i.e., all rows (or, equivalently,
346 columns) must be linearly independent; if either is not true, use
347 `lstsq` for the least-squares best "solution" of the
348 system/equation.
349
350 References
351 ----------
352 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
353 FL, Academic Press, Inc., 1980, pg. 22.
354
355 Examples
356 --------
357 Solve the system of equations ``x0 + 2 * x1 = 1`` and ``3 * x0 + 5 * x1 = 2``:
358
359 >>> a = np.array([[1, 2], [3, 5]])
360 >>> b = np.array([1, 2])
361 >>> x = np.linalg.solve(a, b)
362 >>> x
363 array([-1., 1.])
364
365 Check that the solution is correct:
366
367 >>> np.allclose(np.dot(a, x), b)
368 True
369
370 """
371 a, _ = _makearray(a)
372 _assert_stacked_2d(a)
373 _assert_stacked_square(a)
374 b, wrap = _makearray(b)
375 t, result_t = _commonType(a, b)
376
377 # We use the b = (..., M,) logic, only if the number of extra dimensions
378 # match exactly
379 if b.ndim == a.ndim - 1:
380 gufunc = _umath_linalg.solve1
381 else:
382 gufunc = _umath_linalg.solve
383
384 signature = 'DD->D' if isComplexType(t) else 'dd->d'
385 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
386 r = gufunc(a, b, signature=signature, extobj=extobj)
387
388 return wrap(r.astype(result_t, copy=False))
389
390
391def _tensorinv_dispatcher(a, ind=None):
392 return (a,)
393
394
395@array_function_dispatch(_tensorinv_dispatcher)
396def tensorinv(a, ind=2):
397 """
398 Compute the 'inverse' of an N-dimensional array.
399
400 The result is an inverse for `a` relative to the tensordot operation
401 ``tensordot(a, b, ind)``, i. e., up to floating-point accuracy,
402 ``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the
403 tensordot operation.
404
405 Parameters
406 ----------
407 a : array_like
408 Tensor to 'invert'. Its shape must be 'square', i. e.,
409 ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.
410 ind : int, optional
411 Number of first indices that are involved in the inverse sum.
412 Must be a positive integer, default is 2.
413
414 Returns
415 -------
416 b : ndarray
417 `a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``.
418
419 Raises
420 ------
421 LinAlgError
422 If `a` is singular or not 'square' (in the above sense).
423
424 See Also
425 --------
426 numpy.tensordot, tensorsolve
427
428 Examples
429 --------
430 >>> a = np.eye(4*6)
431 >>> a.shape = (4, 6, 8, 3)
432 >>> ainv = np.linalg.tensorinv(a, ind=2)
433 >>> ainv.shape
434 (8, 3, 4, 6)
435 >>> b = np.random.randn(4, 6)
436 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
437 True
438
439 >>> a = np.eye(4*6)
440 >>> a.shape = (24, 8, 3)
441 >>> ainv = np.linalg.tensorinv(a, ind=1)
442 >>> ainv.shape
443 (8, 3, 24)
444 >>> b = np.random.randn(24)
445 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
446 True
447
448 """
449 a = asarray(a)
450 oldshape = a.shape
451 prod = 1
452 if ind > 0:
453 invshape = oldshape[ind:] + oldshape[:ind]
454 for k in oldshape[ind:]:
455 prod *= k
456 else:
457 raise ValueError("Invalid ind argument.")
458 a = a.reshape(prod, -1)
459 ia = inv(a)
460 return ia.reshape(*invshape)
461
462
463# Matrix inversion
464
465def _unary_dispatcher(a):
466 return (a,)
467
468
469@array_function_dispatch(_unary_dispatcher)
470def inv(a):
471 """
472 Compute the (multiplicative) inverse of a matrix.
473
474 Given a square matrix `a`, return the matrix `ainv` satisfying
475 ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``.
476
477 Parameters
478 ----------
479 a : (..., M, M) array_like
480 Matrix to be inverted.
481
482 Returns
483 -------
484 ainv : (..., M, M) ndarray or matrix
485 (Multiplicative) inverse of the matrix `a`.
486
487 Raises
488 ------
489 LinAlgError
490 If `a` is not square or inversion fails.
491
492 See Also
493 --------
494 scipy.linalg.inv : Similar function in SciPy.
495
496 Notes
497 -----
498
499 .. versionadded:: 1.8.0
500
501 Broadcasting rules apply, see the `numpy.linalg` documentation for
502 details.
503
504 Examples
505 --------
506 >>> from numpy.linalg import inv
507 >>> a = np.array([[1., 2.], [3., 4.]])
508 >>> ainv = inv(a)
509 >>> np.allclose(np.dot(a, ainv), np.eye(2))
510 True
511 >>> np.allclose(np.dot(ainv, a), np.eye(2))
512 True
513
514 If a is a matrix object, then the return value is a matrix as well:
515
516 >>> ainv = inv(np.matrix(a))
517 >>> ainv
518 matrix([[-2. , 1. ],
519 [ 1.5, -0.5]])
520
521 Inverses of several matrices can be computed at once:
522
523 >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
524 >>> inv(a)
525 array([[[-2. , 1. ],
526 [ 1.5 , -0.5 ]],
527 [[-1.25, 0.75],
528 [ 0.75, -0.25]]])
529
530 """
531 a, wrap = _makearray(a)
532 _assert_stacked_2d(a)
533 _assert_stacked_square(a)
534 t, result_t = _commonType(a)
535
536 signature = 'D->D' if isComplexType(t) else 'd->d'
537 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
538 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
539 return wrap(ainv.astype(result_t, copy=False))
540
541
542def _matrix_power_dispatcher(a, n):
543 return (a,)
544
545
546@array_function_dispatch(_matrix_power_dispatcher)
547def matrix_power(a, n):
548 """
549 Raise a square matrix to the (integer) power `n`.
550
551 For positive integers `n`, the power is computed by repeated matrix
552 squarings and matrix multiplications. If ``n == 0``, the identity matrix
553 of the same shape as M is returned. If ``n < 0``, the inverse
554 is computed and then raised to the ``abs(n)``.
555
556 .. note:: Stacks of object matrices are not currently supported.
557
558 Parameters
559 ----------
560 a : (..., M, M) array_like
561 Matrix to be "powered".
562 n : int
563 The exponent can be any integer or long integer, positive,
564 negative, or zero.
565
566 Returns
567 -------
568 a**n : (..., M, M) ndarray or matrix object
569 The return value is the same shape and type as `M`;
570 if the exponent is positive or zero then the type of the
571 elements is the same as those of `M`. If the exponent is
572 negative the elements are floating-point.
573
574 Raises
575 ------
576 LinAlgError
577 For matrices that are not square or that (for negative powers) cannot
578 be inverted numerically.
579
580 Examples
581 --------
582 >>> from numpy.linalg import matrix_power
583 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
584 >>> matrix_power(i, 3) # should = -i
585 array([[ 0, -1],
586 [ 1, 0]])
587 >>> matrix_power(i, 0)
588 array([[1, 0],
589 [0, 1]])
590 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
591 array([[ 0., 1.],
592 [-1., 0.]])
593
594 Somewhat more sophisticated example
595
596 >>> q = np.zeros((4, 4))
597 >>> q[0:2, 0:2] = -i
598 >>> q[2:4, 2:4] = i
599 >>> q # one of the three quaternion units not equal to 1
600 array([[ 0., -1., 0., 0.],
601 [ 1., 0., 0., 0.],
602 [ 0., 0., 0., 1.],
603 [ 0., 0., -1., 0.]])
604 >>> matrix_power(q, 2) # = -np.eye(4)
605 array([[-1., 0., 0., 0.],
606 [ 0., -1., 0., 0.],
607 [ 0., 0., -1., 0.],
608 [ 0., 0., 0., -1.]])
609
610 """
611 a = asanyarray(a)
612 _assert_stacked_2d(a)
613 _assert_stacked_square(a)
614
615 try:
616 n = operator.index(n)
617 except TypeError as e:
618 raise TypeError("exponent must be an integer") from e
619
620 # Fall back on dot for object arrays. Object arrays are not supported by
621 # the current implementation of matmul using einsum
622 if a.dtype != object:
623 fmatmul = matmul
624 elif a.ndim == 2:
625 fmatmul = dot
626 else:
627 raise NotImplementedError(
628 "matrix_power not supported for stacks of object arrays")
629
630 if n == 0:
631 a = empty_like(a)
632 a[...] = eye(a.shape[-2], dtype=a.dtype)
633 return a
634
635 elif n < 0:
636 a = inv(a)
637 n = abs(n)
638
639 # short-cuts.
640 if n == 1:
641 return a
642
643 elif n == 2:
644 return fmatmul(a, a)
645
646 elif n == 3:
647 return fmatmul(fmatmul(a, a), a)
648
649 # Use binary decomposition to reduce the number of matrix multiplications.
650 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
651 # increasing powers of 2, and multiply into the result as needed.
652 z = result = None
653 while n > 0:
654 z = a if z is None else fmatmul(z, z)
655 n, bit = divmod(n, 2)
656 if bit:
657 result = z if result is None else fmatmul(result, z)
658
659 return result
660
661
662# Cholesky decomposition
663
664
665@array_function_dispatch(_unary_dispatcher)
666def cholesky(a):
667 """
668 Cholesky decomposition.
669
670 Return the Cholesky decomposition, `L * L.H`, of the square matrix `a`,
671 where `L` is lower-triangular and .H is the conjugate transpose operator
672 (which is the ordinary transpose if `a` is real-valued). `a` must be
673 Hermitian (symmetric if real-valued) and positive-definite. No
674 checking is performed to verify whether `a` is Hermitian or not.
675 In addition, only the lower-triangular and diagonal elements of `a`
676 are used. Only `L` is actually returned.
677
678 Parameters
679 ----------
680 a : (..., M, M) array_like
681 Hermitian (symmetric if all elements are real), positive-definite
682 input matrix.
683
684 Returns
685 -------
686 L : (..., M, M) array_like
687 Lower-triangular Cholesky factor of `a`. Returns a matrix object if
688 `a` is a matrix object.
689
690 Raises
691 ------
692 LinAlgError
693 If the decomposition fails, for example, if `a` is not
694 positive-definite.
695
696 See Also
697 --------
698 scipy.linalg.cholesky : Similar function in SciPy.
699 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
700 positive-definite matrix.
701 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
702 `scipy.linalg.cho_solve`.
703
704 Notes
705 -----
706
707 .. versionadded:: 1.8.0
708
709 Broadcasting rules apply, see the `numpy.linalg` documentation for
710 details.
711
712 The Cholesky decomposition is often used as a fast way of solving
713
714 .. math:: A \\mathbf{x} = \\mathbf{b}
715
716 (when `A` is both Hermitian/symmetric and positive-definite).
717
718 First, we solve for :math:`\\mathbf{y}` in
719
720 .. math:: L \\mathbf{y} = \\mathbf{b},
721
722 and then for :math:`\\mathbf{x}` in
723
724 .. math:: L.H \\mathbf{x} = \\mathbf{y}.
725
726 Examples
727 --------
728 >>> A = np.array([[1,-2j],[2j,5]])
729 >>> A
730 array([[ 1.+0.j, -0.-2.j],
731 [ 0.+2.j, 5.+0.j]])
732 >>> L = np.linalg.cholesky(A)
733 >>> L
734 array([[1.+0.j, 0.+0.j],
735 [0.+2.j, 1.+0.j]])
736 >>> np.dot(L, L.T.conj()) # verify that L * L.H = A
737 array([[1.+0.j, 0.-2.j],
738 [0.+2.j, 5.+0.j]])
739 >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
740 >>> np.linalg.cholesky(A) # an ndarray object is returned
741 array([[1.+0.j, 0.+0.j],
742 [0.+2.j, 1.+0.j]])
743 >>> # But a matrix object is returned if A is a matrix object
744 >>> np.linalg.cholesky(np.matrix(A))
745 matrix([[ 1.+0.j, 0.+0.j],
746 [ 0.+2.j, 1.+0.j]])
747
748 """
749 extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef)
750 gufunc = _umath_linalg.cholesky_lo
751 a, wrap = _makearray(a)
752 _assert_stacked_2d(a)
753 _assert_stacked_square(a)
754 t, result_t = _commonType(a)
755 signature = 'D->D' if isComplexType(t) else 'd->d'
756 r = gufunc(a, signature=signature, extobj=extobj)
757 return wrap(r.astype(result_t, copy=False))
758
759
760# QR decomposition
761
762def _qr_dispatcher(a, mode=None):
763 return (a,)
764
765
766@array_function_dispatch(_qr_dispatcher)
767def qr(a, mode='reduced'):
768 """
769 Compute the qr factorization of a matrix.
770
771 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is
772 upper-triangular.
773
774 Parameters
775 ----------
776 a : array_like, shape (..., M, N)
777 An array-like object with the dimensionality of at least 2.
778 mode : {'reduced', 'complete', 'r', 'raw'}, optional
779 If K = min(M, N), then
780
781 * 'reduced' : returns q, r with dimensions
782 (..., M, K), (..., K, N) (default)
783 * 'complete' : returns q, r with dimensions (..., M, M), (..., M, N)
784 * 'r' : returns r only with dimensions (..., K, N)
785 * 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)
786
787 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
788 see the notes for more information. The default is 'reduced', and to
789 maintain backward compatibility with earlier versions of numpy both
790 it and the old default 'full' can be omitted. Note that array h
791 returned in 'raw' mode is transposed for calling Fortran. The
792 'economic' mode is deprecated. The modes 'full' and 'economic' may
793 be passed using only the first letter for backwards compatibility,
794 but all others must be spelled out. See the Notes for more
795 explanation.
796
797
798 Returns
799 -------
800 q : ndarray of float or complex, optional
801 A matrix with orthonormal columns. When mode = 'complete' the
802 result is an orthogonal/unitary matrix depending on whether or not
803 a is real/complex. The determinant may be either +/- 1 in that
804 case. In case the number of dimensions in the input array is
805 greater than 2 then a stack of the matrices with above properties
806 is returned.
807 r : ndarray of float or complex, optional
808 The upper-triangular matrix or a stack of upper-triangular
809 matrices if the number of dimensions in the input array is greater
810 than 2.
811 (h, tau) : ndarrays of np.double or np.cdouble, optional
812 The array h contains the Householder reflectors that generate q
813 along with r. The tau array contains scaling factors for the
814 reflectors. In the deprecated 'economic' mode only h is returned.
815
816 Raises
817 ------
818 LinAlgError
819 If factoring fails.
820
821 See Also
822 --------
823 scipy.linalg.qr : Similar function in SciPy.
824 scipy.linalg.rq : Compute RQ decomposition of a matrix.
825
826 Notes
827 -----
828 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
829 ``dorgqr``, and ``zungqr``.
830
831 For more information on the qr factorization, see for example:
832 https://en.wikipedia.org/wiki/QR_factorization
833
834 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
835 `a` is of type `matrix`, all the return values will be matrices too.
836
837 New 'reduced', 'complete', and 'raw' options for mode were added in
838 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In
839 addition the options 'full' and 'economic' were deprecated. Because
840 'full' was the previous default and 'reduced' is the new default,
841 backward compatibility can be maintained by letting `mode` default.
842 The 'raw' option was added so that LAPACK routines that can multiply
843 arrays by q using the Householder reflectors can be used. Note that in
844 this case the returned arrays are of type np.double or np.cdouble and
845 the h array is transposed to be FORTRAN compatible. No routines using
846 the 'raw' return are currently exposed by numpy, but some are available
847 in lapack_lite and just await the necessary work.
848
849 Examples
850 --------
851 >>> a = np.random.randn(9, 6)
852 >>> q, r = np.linalg.qr(a)
853 >>> np.allclose(a, np.dot(q, r)) # a does equal qr
854 True
855 >>> r2 = np.linalg.qr(a, mode='r')
856 >>> np.allclose(r, r2) # mode='r' returns the same r as mode='full'
857 True
858 >>> a = np.random.normal(size=(3, 2, 2)) # Stack of 2 x 2 matrices as input
859 >>> q, r = np.linalg.qr(a)
860 >>> q.shape
861 (3, 2, 2)
862 >>> r.shape
863 (3, 2, 2)
864 >>> np.allclose(a, np.matmul(q, r))
865 True
866
867 Example illustrating a common use of `qr`: solving of least squares
868 problems
869
870 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for
871 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
872 and you'll see that it should be y0 = 0, m = 1.) The answer is provided
873 by solving the over-determined matrix equation ``Ax = b``, where::
874
875 A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
876 x = array([[y0], [m]])
877 b = array([[1], [0], [2], [1]])
878
879 If A = qr such that q is orthonormal (which is always possible via
880 Gram-Schmidt), then ``x = inv(r) * (q.T) * b``. (In numpy practice,
881 however, we simply use `lstsq`.)
882
883 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
884 >>> A
885 array([[0, 1],
886 [1, 1],
887 [1, 1],
888 [2, 1]])
889 >>> b = np.array([1, 2, 2, 3])
890 >>> q, r = np.linalg.qr(A)
891 >>> p = np.dot(q.T, b)
892 >>> np.dot(np.linalg.inv(r), p)
893 array([ 1., 1.])
894
895 """
896 if mode not in ('reduced', 'complete', 'r', 'raw'):
897 if mode in ('f', 'full'):
898 # 2013-04-01, 1.8
899 msg = "".join((
900 "The 'full' option is deprecated in favor of 'reduced'.\n",
901 "For backward compatibility let mode default."))
902 warnings.warn(msg, DeprecationWarning, stacklevel=3)
903 mode = 'reduced'
904 elif mode in ('e', 'economic'):
905 # 2013-04-01, 1.8
906 msg = "The 'economic' option is deprecated."
907 warnings.warn(msg, DeprecationWarning, stacklevel=3)
908 mode = 'economic'
909 else:
910 raise ValueError(f"Unrecognized mode '{mode}'")
911
912 a, wrap = _makearray(a)
913 _assert_stacked_2d(a)
914 m, n = a.shape[-2:]
915 t, result_t = _commonType(a)
916 a = a.astype(t, copy=True)
917 a = _to_native_byte_order(a)
918 mn = min(m, n)
919
920 if m <= n:
921 gufunc = _umath_linalg.qr_r_raw_m
922 else:
923 gufunc = _umath_linalg.qr_r_raw_n
924
925 signature = 'D->D' if isComplexType(t) else 'd->d'
926 extobj = get_linalg_error_extobj(_raise_linalgerror_qr)
927 tau = gufunc(a, signature=signature, extobj=extobj)
928
929 # handle modes that don't return q
930 if mode == 'r':
931 r = triu(a[..., :mn, :])
932 r = r.astype(result_t, copy=False)
933 return wrap(r)
934
935 if mode == 'raw':
936 q = transpose(a)
937 q = q.astype(result_t, copy=False)
938 tau = tau.astype(result_t, copy=False)
939 return wrap(q), tau
940
941 if mode == 'economic':
942 a = a.astype(result_t, copy=False)
943 return wrap(a)
944
945 # mc is the number of columns in the resulting q
946 # matrix. If the mode is complete then it is
947 # same as number of rows, and if the mode is reduced,
948 # then it is the minimum of number of rows and columns.
949 if mode == 'complete' and m > n:
950 mc = m
951 gufunc = _umath_linalg.qr_complete
952 else:
953 mc = mn
954 gufunc = _umath_linalg.qr_reduced
955
956 signature = 'DD->D' if isComplexType(t) else 'dd->d'
957 extobj = get_linalg_error_extobj(_raise_linalgerror_qr)
958 q = gufunc(a, tau, signature=signature, extobj=extobj)
959 r = triu(a[..., :mc, :])
960
961 q = q.astype(result_t, copy=False)
962 r = r.astype(result_t, copy=False)
963
964 return wrap(q), wrap(r)
965
966# Eigenvalues
967
968
969@array_function_dispatch(_unary_dispatcher)
970def eigvals(a):
971 """
972 Compute the eigenvalues of a general matrix.
973
974 Main difference between `eigvals` and `eig`: the eigenvectors aren't
975 returned.
976
977 Parameters
978 ----------
979 a : (..., M, M) array_like
980 A complex- or real-valued matrix whose eigenvalues will be computed.
981
982 Returns
983 -------
984 w : (..., M,) ndarray
985 The eigenvalues, each repeated according to its multiplicity.
986 They are not necessarily ordered, nor are they necessarily
987 real for real matrices.
988
989 Raises
990 ------
991 LinAlgError
992 If the eigenvalue computation does not converge.
993
994 See Also
995 --------
996 eig : eigenvalues and right eigenvectors of general arrays
997 eigvalsh : eigenvalues of real symmetric or complex Hermitian
998 (conjugate symmetric) arrays.
999 eigh : eigenvalues and eigenvectors of real symmetric or complex
1000 Hermitian (conjugate symmetric) arrays.
1001 scipy.linalg.eigvals : Similar function in SciPy.
1002
1003 Notes
1004 -----
1005
1006 .. versionadded:: 1.8.0
1007
1008 Broadcasting rules apply, see the `numpy.linalg` documentation for
1009 details.
1010
1011 This is implemented using the ``_geev`` LAPACK routines which compute
1012 the eigenvalues and eigenvectors of general square arrays.
1013
1014 Examples
1015 --------
1016 Illustration, using the fact that the eigenvalues of a diagonal matrix
1017 are its diagonal elements, that multiplying a matrix on the left
1018 by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
1019 of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
1020 if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
1021 ``A``:
1022
1023 >>> from numpy import linalg as LA
1024 >>> x = np.random.random()
1025 >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
1026 >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :])
1027 (1.0, 1.0, 0.0)
1028
1029 Now multiply a diagonal matrix by ``Q`` on one side and by ``Q.T`` on the other:
1030
1031 >>> D = np.diag((-1,1))
1032 >>> LA.eigvals(D)
1033 array([-1., 1.])
1034 >>> A = np.dot(Q, D)
1035 >>> A = np.dot(A, Q.T)
1036 >>> LA.eigvals(A)
1037 array([ 1., -1.]) # random
1038
1039 """
1040 a, wrap = _makearray(a)
1041 _assert_stacked_2d(a)
1042 _assert_stacked_square(a)
1043 _assert_finite(a)
1044 t, result_t = _commonType(a)
1045
1046 extobj = get_linalg_error_extobj(
1047 _raise_linalgerror_eigenvalues_nonconvergence)
1048 signature = 'D->D' if isComplexType(t) else 'd->D'
1049 w = _umath_linalg.eigvals(a, signature=signature, extobj=extobj)
1050
1051 if not isComplexType(t):
1052 if all(w.imag == 0):
1053 w = w.real
1054 result_t = _realType(result_t)
1055 else:
1056 result_t = _complexType(result_t)
1057
1058 return w.astype(result_t, copy=False)
1059
1060
1061def _eigvalsh_dispatcher(a, UPLO=None):
1062 return (a,)
1063
1064
1065@array_function_dispatch(_eigvalsh_dispatcher)
1066def eigvalsh(a, UPLO='L'):
1067 """
1068 Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
1069
1070 Main difference from eigh: the eigenvectors are not computed.
1071
1072 Parameters
1073 ----------
1074 a : (..., M, M) array_like
1075 A complex- or real-valued matrix whose eigenvalues are to be
1076 computed.
1077 UPLO : {'L', 'U'}, optional
1078 Specifies whether the calculation is done with the lower triangular
1079 part of `a` ('L', default) or the upper triangular part ('U').
1080 Irrespective of this value only the real parts of the diagonal will
1081 be considered in the computation to preserve the notion of a Hermitian
1082 matrix. It therefore follows that the imaginary part of the diagonal
1083 will always be treated as zero.
1084
1085 Returns
1086 -------
1087 w : (..., M,) ndarray
1088 The eigenvalues in ascending order, each repeated according to
1089 its multiplicity.
1090
1091 Raises
1092 ------
1093 LinAlgError
1094 If the eigenvalue computation does not converge.
1095
1096 See Also
1097 --------
1098 eigh : eigenvalues and eigenvectors of real symmetric or complex Hermitian
1099 (conjugate symmetric) arrays.
1100 eigvals : eigenvalues of general real or complex arrays.
1101 eig : eigenvalues and right eigenvectors of general real or complex
1102 arrays.
1103 scipy.linalg.eigvalsh : Similar function in SciPy.
1104
1105 Notes
1106 -----
1107
1108 .. versionadded:: 1.8.0
1109
1110 Broadcasting rules apply, see the `numpy.linalg` documentation for
1111 details.
1112
1113 The eigenvalues are computed using LAPACK routines ``_syevd``, ``_heevd``.
1114
1115 Examples
1116 --------
1117 >>> from numpy import linalg as LA
1118 >>> a = np.array([[1, -2j], [2j, 5]])
1119 >>> LA.eigvalsh(a)
1120 array([ 0.17157288, 5.82842712]) # may vary
1121
1122 >>> # demonstrate the treatment of the imaginary part of the diagonal
1123 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1124 >>> a
1125 array([[5.+2.j, 9.-2.j],
1126 [0.+2.j, 2.-1.j]])
1127 >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals()
1128 >>> # with:
1129 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1130 >>> b
1131 array([[5.+0.j, 0.-2.j],
1132 [0.+2.j, 2.+0.j]])
1133 >>> wa = LA.eigvalsh(a)
1134 >>> wb = LA.eigvals(b)
1135 >>> wa; wb
1136 array([1., 6.])
1137 array([6.+0.j, 1.+0.j])
1138
1139 """
1140 UPLO = UPLO.upper()
1141 if UPLO not in ('L', 'U'):
1142 raise ValueError("UPLO argument must be 'L' or 'U'")
1143
1144 extobj = get_linalg_error_extobj(
1145 _raise_linalgerror_eigenvalues_nonconvergence)
1146 if UPLO == 'L':
1147 gufunc = _umath_linalg.eigvalsh_lo
1148 else:
1149 gufunc = _umath_linalg.eigvalsh_up
1150
1151 a, wrap = _makearray(a)
1152 _assert_stacked_2d(a)
1153 _assert_stacked_square(a)
1154 t, result_t = _commonType(a)
1155 signature = 'D->d' if isComplexType(t) else 'd->d'
1156 w = gufunc(a, signature=signature, extobj=extobj)
1157 return w.astype(_realType(result_t), copy=False)
1158
1159def _convertarray(a):
1160 t, result_t = _commonType(a)
1161 a = a.astype(t).T.copy()
1162 return a, t, result_t
1163
1164
1165# Eigenvectors
1166
1167
1168@array_function_dispatch(_unary_dispatcher)
1169def eig(a):
1170 """
1171 Compute the eigenvalues and right eigenvectors of a square array.
1172
1173 Parameters
1174 ----------
1175 a : (..., M, M) array
1176 Matrices for which the eigenvalues and right eigenvectors will
1177 be computed
1178
1179 Returns
1180 -------
1181 w : (..., M) array
1182 The eigenvalues, each repeated according to its multiplicity.
1183 The eigenvalues are not necessarily ordered. The resulting
1184 array will be of complex type, unless the imaginary part is
1185 zero in which case it will be cast to a real type. When `a`
1186 is real the resulting eigenvalues will be real (0 imaginary
1187 part) or occur in conjugate pairs
1188
1189 v : (..., M, M) array
1190 The normalized (unit "length") eigenvectors, such that the
1191 column ``v[:,i]`` is the eigenvector corresponding to the
1192 eigenvalue ``w[i]``.
1193
1194 Raises
1195 ------
1196 LinAlgError
1197 If the eigenvalue computation does not converge.
1198
1199 See Also
1200 --------
1201 eigvals : eigenvalues of a non-symmetric array.
1202 eigh : eigenvalues and eigenvectors of a real symmetric or complex
1203 Hermitian (conjugate symmetric) array.
1204 eigvalsh : eigenvalues of a real symmetric or complex Hermitian
1205 (conjugate symmetric) array.
1206 scipy.linalg.eig : Similar function in SciPy that also solves the
1207 generalized eigenvalue problem.
1208 scipy.linalg.schur : Best choice for unitary and other non-Hermitian
1209 normal matrices.
1210
1211 Notes
1212 -----
1213
1214 .. versionadded:: 1.8.0
1215
1216 Broadcasting rules apply, see the `numpy.linalg` documentation for
1217 details.
1218
1219 This is implemented using the ``_geev`` LAPACK routines which compute
1220 the eigenvalues and eigenvectors of general square arrays.
1221
1222 The number `w` is an eigenvalue of `a` if there exists a vector
1223 `v` such that ``a @ v = w * v``. Thus, the arrays `a`, `w`, and
1224 `v` satisfy the equations ``a @ v[:,i] = w[i] * v[:,i]``
1225 for :math:`i \\in \\{0,...,M-1\\}`.
1226
1227 The array `v` of eigenvectors may not be of maximum rank, that is, some
1228 of the columns may be linearly dependent, although round-off error may
1229 obscure that fact. If the eigenvalues are all different, then theoretically
1230 the eigenvectors are linearly independent and `a` can be diagonalized by
1231 a similarity transformation using `v`, i.e, ``inv(v) @ a @ v`` is diagonal.
1232
1233 For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
1234 is preferred because the matrix `v` is guaranteed to be unitary, which is
1235 not the case when using `eig`. The Schur factorization produces an
1236 upper triangular matrix rather than a diagonal matrix, but for normal
1237 matrices only the diagonal of the upper triangular matrix is needed, the
1238 rest is roundoff error.
1239
1240 Finally, it is emphasized that `v` consists of the *right* (as in
1241 right-hand side) eigenvectors of `a`. A vector `y` satisfying
1242 ``y.T @ a = z * y.T`` for some number `z` is called a *left*
1243 eigenvector of `a`, and, in general, the left and right eigenvectors
1244 of a matrix are not necessarily the (perhaps conjugate) transposes
1245 of each other.
1246
1247 References
1248 ----------
1249 G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
1250 Academic Press, Inc., 1980, Various pp.
1251
1252 Examples
1253 --------
1254 >>> from numpy import linalg as LA
1255
1256 (Almost) trivial example with real e-values and e-vectors.
1257
1258 >>> w, v = LA.eig(np.diag((1, 2, 3)))
1259 >>> w; v
1260 array([1., 2., 3.])
1261 array([[1., 0., 0.],
1262 [0., 1., 0.],
1263 [0., 0., 1.]])
1264
1265 Real matrix possessing complex e-values and e-vectors; note that the
1266 e-values are complex conjugates of each other.
1267
1268 >>> w, v = LA.eig(np.array([[1, -1], [1, 1]]))
1269 >>> w; v
1270 array([1.+1.j, 1.-1.j])
1271 array([[0.70710678+0.j , 0.70710678-0.j ],
1272 [0. -0.70710678j, 0. +0.70710678j]])
1273
1274 Complex-valued matrix with real e-values (but complex-valued e-vectors);
1275 note that ``a.conj().T == a``, i.e., `a` is Hermitian.
1276
1277 >>> a = np.array([[1, 1j], [-1j, 1]])
1278 >>> w, v = LA.eig(a)
1279 >>> w; v
1280 array([2.+0.j, 0.+0.j])
1281 array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary
1282 [ 0.70710678+0.j , -0. +0.70710678j]])
1283
1284 Be careful about round-off error!
1285
1286 >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]])
1287 >>> # Theor. e-values are 1 +/- 1e-9
1288 >>> w, v = LA.eig(a)
1289 >>> w; v
1290 array([1., 1.])
1291 array([[1., 0.],
1292 [0., 1.]])
1293
1294 """
1295 a, wrap = _makearray(a)
1296 _assert_stacked_2d(a)
1297 _assert_stacked_square(a)
1298 _assert_finite(a)
1299 t, result_t = _commonType(a)
1300
1301 extobj = get_linalg_error_extobj(
1302 _raise_linalgerror_eigenvalues_nonconvergence)
1303 signature = 'D->DD' if isComplexType(t) else 'd->DD'
1304 w, vt = _umath_linalg.eig(a, signature=signature, extobj=extobj)
1305
1306 if not isComplexType(t) and all(w.imag == 0.0):
1307 w = w.real
1308 vt = vt.real
1309 result_t = _realType(result_t)
1310 else:
1311 result_t = _complexType(result_t)
1312
1313 vt = vt.astype(result_t, copy=False)
1314 return w.astype(result_t, copy=False), wrap(vt)
1315
1316
1317@array_function_dispatch(_eigvalsh_dispatcher)
1318def eigh(a, UPLO='L'):
1319 """
1320 Return the eigenvalues and eigenvectors of a complex Hermitian
1321 (conjugate symmetric) or a real symmetric matrix.
1322
1323 Returns two objects, a 1-D array containing the eigenvalues of `a`, and
1324 a 2-D square array or matrix (depending on the input type) of the
1325 corresponding eigenvectors (in columns).
1326
1327 Parameters
1328 ----------
1329 a : (..., M, M) array
1330 Hermitian or real symmetric matrices whose eigenvalues and
1331 eigenvectors are to be computed.
1332 UPLO : {'L', 'U'}, optional
1333 Specifies whether the calculation is done with the lower triangular
1334 part of `a` ('L', default) or the upper triangular part ('U').
1335 Irrespective of this value only the real parts of the diagonal will
1336 be considered in the computation to preserve the notion of a Hermitian
1337 matrix. It therefore follows that the imaginary part of the diagonal
1338 will always be treated as zero.
1339
1340 Returns
1341 -------
1342 w : (..., M) ndarray
1343 The eigenvalues in ascending order, each repeated according to
1344 its multiplicity.
1345 v : {(..., M, M) ndarray, (..., M, M) matrix}
1346 The column ``v[:, i]`` is the normalized eigenvector corresponding
1347 to the eigenvalue ``w[i]``. Will return a matrix object if `a` is
1348 a matrix object.
1349
1350 Raises
1351 ------
1352 LinAlgError
1353 If the eigenvalue computation does not converge.
1354
1355 See Also
1356 --------
1357 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1358 (conjugate symmetric) arrays.
1359 eig : eigenvalues and right eigenvectors for non-symmetric arrays.
1360 eigvals : eigenvalues of non-symmetric arrays.
1361 scipy.linalg.eigh : Similar function in SciPy (but also solves the
1362 generalized eigenvalue problem).
1363
1364 Notes
1365 -----
1366
1367 .. versionadded:: 1.8.0
1368
1369 Broadcasting rules apply, see the `numpy.linalg` documentation for
1370 details.
1371
1372 The eigenvalues/eigenvectors are computed using LAPACK routines ``_syevd``,
1373 ``_heevd``.
1374
1375 The eigenvalues of real symmetric or complex Hermitian matrices are
1376 always real. [1]_ The array `v` of (column) eigenvectors is unitary
1377 and `a`, `w`, and `v` satisfy the equations
1378 ``dot(a, v[:, i]) = w[i] * v[:, i]``.
1379
1380 References
1381 ----------
1382 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1383 FL, Academic Press, Inc., 1980, pg. 222.
1384
1385 Examples
1386 --------
1387 >>> from numpy import linalg as LA
1388 >>> a = np.array([[1, -2j], [2j, 5]])
1389 >>> a
1390 array([[ 1.+0.j, -0.-2.j],
1391 [ 0.+2.j, 5.+0.j]])
1392 >>> w, v = LA.eigh(a)
1393 >>> w; v
1394 array([0.17157288, 5.82842712])
1395 array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1396 [ 0. +0.38268343j, 0. -0.92387953j]])
1397
1398 >>> np.dot(a, v[:, 0]) - w[0] * v[:, 0] # verify 1st e-val/vec pair
1399 array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j])
1400 >>> np.dot(a, v[:, 1]) - w[1] * v[:, 1] # verify 2nd e-val/vec pair
1401 array([0.+0.j, 0.+0.j])
1402
1403 >>> A = np.matrix(a) # what happens if input is a matrix object
1404 >>> A
1405 matrix([[ 1.+0.j, -0.-2.j],
1406 [ 0.+2.j, 5.+0.j]])
1407 >>> w, v = LA.eigh(A)
1408 >>> w; v
1409 array([0.17157288, 5.82842712])
1410 matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1411 [ 0. +0.38268343j, 0. -0.92387953j]])
1412
1413 >>> # demonstrate the treatment of the imaginary part of the diagonal
1414 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1415 >>> a
1416 array([[5.+2.j, 9.-2.j],
1417 [0.+2.j, 2.-1.j]])
1418 >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with:
1419 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1420 >>> b
1421 array([[5.+0.j, 0.-2.j],
1422 [0.+2.j, 2.+0.j]])
1423 >>> wa, va = LA.eigh(a)
1424 >>> wb, vb = LA.eig(b)
1425 >>> wa; wb
1426 array([1., 6.])
1427 array([6.+0.j, 1.+0.j])
1428 >>> va; vb
1429 array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary
1430 [ 0. +0.89442719j, 0. -0.4472136j ]])
1431 array([[ 0.89442719+0.j , -0. +0.4472136j],
1432 [-0. +0.4472136j, 0.89442719+0.j ]])
1433 """
1434 UPLO = UPLO.upper()
1435 if UPLO not in ('L', 'U'):
1436 raise ValueError("UPLO argument must be 'L' or 'U'")
1437
1438 a, wrap = _makearray(a)
1439 _assert_stacked_2d(a)
1440 _assert_stacked_square(a)
1441 t, result_t = _commonType(a)
1442
1443 extobj = get_linalg_error_extobj(
1444 _raise_linalgerror_eigenvalues_nonconvergence)
1445 if UPLO == 'L':
1446 gufunc = _umath_linalg.eigh_lo
1447 else:
1448 gufunc = _umath_linalg.eigh_up
1449
1450 signature = 'D->dD' if isComplexType(t) else 'd->dd'
1451 w, vt = gufunc(a, signature=signature, extobj=extobj)
1452 w = w.astype(_realType(result_t), copy=False)
1453 vt = vt.astype(result_t, copy=False)
1454 return w, wrap(vt)
1455
1456
1457# Singular value decomposition
1458
1459def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None):
1460 return (a,)
1461
1462
1463@array_function_dispatch(_svd_dispatcher)
1464def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
1465 """
1466 Singular Value Decomposition.
1467
1468 When `a` is a 2D array, and ``full_matrices=False``, then it is
1469 factorized as ``u @ np.diag(s) @ vh = (u * s) @ vh``, where
1470 `u` and the Hermitian transpose of `vh` are 2D arrays with
1471 orthonormal columns and `s` is a 1D array of `a`'s singular
1472 values. When `a` is higher-dimensional, SVD is applied in
1473 stacked mode as explained below.
1474
1475 Parameters
1476 ----------
1477 a : (..., M, N) array_like
1478 A real or complex array with ``a.ndim >= 2``.
1479 full_matrices : bool, optional
1480 If True (default), `u` and `vh` have the shapes ``(..., M, M)`` and
1481 ``(..., N, N)``, respectively. Otherwise, the shapes are
1482 ``(..., M, K)`` and ``(..., K, N)``, respectively, where
1483 ``K = min(M, N)``.
1484 compute_uv : bool, optional
1485 Whether or not to compute `u` and `vh` in addition to `s`. True
1486 by default.
1487 hermitian : bool, optional
1488 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1489 enabling a more efficient method for finding singular values.
1490 Defaults to False.
1491
1492 .. versionadded:: 1.17.0
1493
1494 Returns
1495 -------
1496 u : { (..., M, M), (..., M, K) } array
1497 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1498 size as those of the input `a`. The size of the last two dimensions
1499 depends on the value of `full_matrices`. Only returned when
1500 `compute_uv` is True.
1501 s : (..., K) array
1502 Vector(s) with the singular values, within each vector sorted in
1503 descending order. The first ``a.ndim - 2`` dimensions have the same
1504 size as those of the input `a`.
1505 vh : { (..., N, N), (..., K, N) } array
1506 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1507 size as those of the input `a`. The size of the last two dimensions
1508 depends on the value of `full_matrices`. Only returned when
1509 `compute_uv` is True.
1510
1511 Raises
1512 ------
1513 LinAlgError
1514 If SVD computation does not converge.
1515
1516 See Also
1517 --------
1518 scipy.linalg.svd : Similar function in SciPy.
1519 scipy.linalg.svdvals : Compute singular values of a matrix.
1520
1521 Notes
1522 -----
1523
1524 .. versionchanged:: 1.8.0
1525 Broadcasting rules apply, see the `numpy.linalg` documentation for
1526 details.
1527
1528 The decomposition is performed using LAPACK routine ``_gesdd``.
1529
1530 SVD is usually described for the factorization of a 2D matrix :math:`A`.
1531 The higher-dimensional case will be discussed below. In the 2D case, SVD is
1532 written as :math:`A = U S V^H`, where :math:`A = a`, :math:`U= u`,
1533 :math:`S= \\mathtt{np.diag}(s)` and :math:`V^H = vh`. The 1D array `s`
1534 contains the singular values of `a` and `u` and `vh` are unitary. The rows
1535 of `vh` are the eigenvectors of :math:`A^H A` and the columns of `u` are
1536 the eigenvectors of :math:`A A^H`. In both cases the corresponding
1537 (possibly non-zero) eigenvalues are given by ``s**2``.
1538
1539 If `a` has more than two dimensions, then broadcasting rules apply, as
1540 explained in :ref:`routines.linalg-broadcasting`. This means that SVD is
1541 working in "stacked" mode: it iterates over all indices of the first
1542 ``a.ndim - 2`` dimensions and for each combination SVD is applied to the
1543 last two indices. The matrix `a` can be reconstructed from the
1544 decomposition with either ``(u * s[..., None, :]) @ vh`` or
1545 ``u @ (s[..., None] * vh)``. (The ``@`` operator can be replaced by the
1546 function ``np.matmul`` for python versions below 3.5.)
1547
1548 If `a` is a ``matrix`` object (as opposed to an ``ndarray``), then so are
1549 all the return values.
1550
1551 Examples
1552 --------
1553 >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
1554 >>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)
1555
1556 Reconstruction based on full SVD, 2D case:
1557
1558 >>> u, s, vh = np.linalg.svd(a, full_matrices=True)
1559 >>> u.shape, s.shape, vh.shape
1560 ((9, 9), (6,), (6, 6))
1561 >>> np.allclose(a, np.dot(u[:, :6] * s, vh))
1562 True
1563 >>> smat = np.zeros((9, 6), dtype=complex)
1564 >>> smat[:6, :6] = np.diag(s)
1565 >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
1566 True
1567
1568 Reconstruction based on reduced SVD, 2D case:
1569
1570 >>> u, s, vh = np.linalg.svd(a, full_matrices=False)
1571 >>> u.shape, s.shape, vh.shape
1572 ((9, 6), (6,), (6, 6))
1573 >>> np.allclose(a, np.dot(u * s, vh))
1574 True
1575 >>> smat = np.diag(s)
1576 >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
1577 True
1578
1579 Reconstruction based on full SVD, 4D case:
1580
1581 >>> u, s, vh = np.linalg.svd(b, full_matrices=True)
1582 >>> u.shape, s.shape, vh.shape
1583 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
1584 >>> np.allclose(b, np.matmul(u[..., :3] * s[..., None, :], vh))
1585 True
1586 >>> np.allclose(b, np.matmul(u[..., :3], s[..., None] * vh))
1587 True
1588
1589 Reconstruction based on reduced SVD, 4D case:
1590
1591 >>> u, s, vh = np.linalg.svd(b, full_matrices=False)
1592 >>> u.shape, s.shape, vh.shape
1593 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
1594 >>> np.allclose(b, np.matmul(u * s[..., None, :], vh))
1595 True
1596 >>> np.allclose(b, np.matmul(u, s[..., None] * vh))
1597 True
1598
1599 """
1600 import numpy as _nx
1601 a, wrap = _makearray(a)
1602
1603 if hermitian:
1604 # note: lapack svd returns eigenvalues with s ** 2 sorted descending,
1605 # but eig returns s sorted ascending, so we re-order the eigenvalues
1606 # and related arrays to have the correct order
1607 if compute_uv:
1608 s, u = eigh(a)
1609 sgn = sign(s)
1610 s = abs(s)
1611 sidx = argsort(s)[..., ::-1]
1612 sgn = _nx.take_along_axis(sgn, sidx, axis=-1)
1613 s = _nx.take_along_axis(s, sidx, axis=-1)
1614 u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1)
1615 # singular values are unsigned, move the sign into v
1616 vt = transpose(u * sgn[..., None, :]).conjugate()
1617 return wrap(u), s, wrap(vt)
1618 else:
1619 s = eigvalsh(a)
1620 s = abs(s)
1621 return sort(s)[..., ::-1]
1622
1623 _assert_stacked_2d(a)
1624 t, result_t = _commonType(a)
1625
1626 extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence)
1627
1628 m, n = a.shape[-2:]
1629 if compute_uv:
1630 if full_matrices:
1631 if m < n:
1632 gufunc = _umath_linalg.svd_m_f
1633 else:
1634 gufunc = _umath_linalg.svd_n_f
1635 else:
1636 if m < n:
1637 gufunc = _umath_linalg.svd_m_s
1638 else:
1639 gufunc = _umath_linalg.svd_n_s
1640
1641 signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
1642 u, s, vh = gufunc(a, signature=signature, extobj=extobj)
1643 u = u.astype(result_t, copy=False)
1644 s = s.astype(_realType(result_t), copy=False)
1645 vh = vh.astype(result_t, copy=False)
1646 return wrap(u), s, wrap(vh)
1647 else:
1648 if m < n:
1649 gufunc = _umath_linalg.svd_m
1650 else:
1651 gufunc = _umath_linalg.svd_n
1652
1653 signature = 'D->d' if isComplexType(t) else 'd->d'
1654 s = gufunc(a, signature=signature, extobj=extobj)
1655 s = s.astype(_realType(result_t), copy=False)
1656 return s
1657
1658
1659def _cond_dispatcher(x, p=None):
1660 return (x,)
1661
1662
1663@array_function_dispatch(_cond_dispatcher)
1664def cond(x, p=None):
1665 """
1666 Compute the condition number of a matrix.
1667
1668 This function is capable of returning the condition number using
1669 one of seven different norms, depending on the value of `p` (see
1670 Parameters below).
1671
1672 Parameters
1673 ----------
1674 x : (..., M, N) array_like
1675 The matrix whose condition number is sought.
1676 p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
1677 Order of the norm used in the condition number computation:
1678
1679 ===== ============================
1680 p norm for matrices
1681 ===== ============================
1682 None 2-norm, computed directly using the ``SVD``
1683 'fro' Frobenius norm
1684 inf max(sum(abs(x), axis=1))
1685 -inf min(sum(abs(x), axis=1))
1686 1 max(sum(abs(x), axis=0))
1687 -1 min(sum(abs(x), axis=0))
1688 2 2-norm (largest sing. value)
1689 -2 smallest singular value
1690 ===== ============================
1691
1692 inf means the `numpy.inf` object, and the Frobenius norm is
1693 the root-of-sum-of-squares norm.
1694
1695 Returns
1696 -------
1697 c : {float, inf}
1698 The condition number of the matrix. May be infinite.
1699
1700 See Also
1701 --------
1702 numpy.linalg.norm
1703
1704 Notes
1705 -----
1706 The condition number of `x` is defined as the norm of `x` times the
1707 norm of the inverse of `x` [1]_; the norm can be the usual L2-norm
1708 (root-of-sum-of-squares) or one of a number of other matrix norms.
1709
1710 References
1711 ----------
1712 .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL,
1713 Academic Press, Inc., 1980, pg. 285.
1714
1715 Examples
1716 --------
1717 >>> from numpy import linalg as LA
1718 >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
1719 >>> a
1720 array([[ 1, 0, -1],
1721 [ 0, 1, 0],
1722 [ 1, 0, 1]])
1723 >>> LA.cond(a)
1724 1.4142135623730951
1725 >>> LA.cond(a, 'fro')
1726 3.1622776601683795
1727 >>> LA.cond(a, np.inf)
1728 2.0
1729 >>> LA.cond(a, -np.inf)
1730 1.0
1731 >>> LA.cond(a, 1)
1732 2.0
1733 >>> LA.cond(a, -1)
1734 1.0
1735 >>> LA.cond(a, 2)
1736 1.4142135623730951
1737 >>> LA.cond(a, -2)
1738 0.70710678118654746 # may vary
1739 >>> min(LA.svd(a, compute_uv=False))*min(LA.svd(LA.inv(a), compute_uv=False))
1740 0.70710678118654746 # may vary
1741
1742 """
1743 x = asarray(x) # in case we have a matrix
1744 if _is_empty_2d(x):
1745 raise LinAlgError("cond is not defined on empty arrays")
1746 if p is None or p == 2 or p == -2:
1747 s = svd(x, compute_uv=False)
1748 with errstate(all='ignore'):
1749 if p == -2:
1750 r = s[..., -1] / s[..., 0]
1751 else:
1752 r = s[..., 0] / s[..., -1]
1753 else:
1754 # Call inv(x) ignoring errors. The result array will
1755 # contain nans in the entries where inversion failed.
1756 _assert_stacked_2d(x)
1757 _assert_stacked_square(x)
1758 t, result_t = _commonType(x)
1759 signature = 'D->D' if isComplexType(t) else 'd->d'
1760 with errstate(all='ignore'):
1761 invx = _umath_linalg.inv(x, signature=signature)
1762 r = norm(x, p, axis=(-2, -1)) * norm(invx, p, axis=(-2, -1))
1763 r = r.astype(result_t, copy=False)
1764
1765 # Convert nans to infs unless the original array had nan entries
1766 r = asarray(r)
1767 nan_mask = isnan(r)
1768 if nan_mask.any():
1769 nan_mask &= ~isnan(x).any(axis=(-2, -1))
1770 if r.ndim > 0:
1771 r[nan_mask] = Inf
1772 elif nan_mask:
1773 r[()] = Inf
1774
1775 # Convention is to return scalars instead of 0d arrays
1776 if r.ndim == 0:
1777 r = r[()]
1778
1779 return r
1780
1781
1782def _matrix_rank_dispatcher(A, tol=None, hermitian=None):
1783 return (A,)
1784
1785
1786@array_function_dispatch(_matrix_rank_dispatcher)
1787def matrix_rank(A, tol=None, hermitian=False):
1788 """
1789 Return matrix rank of array using SVD method
1790
1791 Rank of the array is the number of singular values of the array that are
1792 greater than `tol`.
1793
1794 .. versionchanged:: 1.14
1795 Can now operate on stacks of matrices
1796
1797 Parameters
1798 ----------
1799 A : {(M,), (..., M, N)} array_like
1800 Input vector or stack of matrices.
1801 tol : (...) array_like, float, optional
1802 Threshold below which SVD values are considered zero. If `tol` is
1803 None, and ``S`` is an array with singular values for `M`, and
1804 ``eps`` is the epsilon value for datatype of ``S``, then `tol` is
1805 set to ``S.max() * max(M, N) * eps``.
1806
1807 .. versionchanged:: 1.14
1808 Broadcasted against the stack of matrices
1809 hermitian : bool, optional
1810 If True, `A` is assumed to be Hermitian (symmetric if real-valued),
1811 enabling a more efficient method for finding singular values.
1812 Defaults to False.
1813
1814 .. versionadded:: 1.14
1815
1816 Returns
1817 -------
1818 rank : (...) array_like
1819 Rank of A.
1820
1821 Notes
1822 -----
1823 The default threshold to detect rank deficiency is a test on the magnitude
1824 of the singular values of `A`. By default, we identify singular values less
1825 than ``S.max() * max(M, N) * eps`` as indicating rank deficiency (with
1826 the symbols defined above). This is the algorithm MATLAB uses [1]. It also
1827 appears in *Numerical recipes* in the discussion of SVD solutions for linear
1828 least squares [2].
1829
1830 This default threshold is designed to detect rank deficiency accounting for
1831 the numerical errors of the SVD computation. Imagine that there is a column
1832 in `A` that is an exact (in floating point) linear combination of other
1833 columns in `A`. Computing the SVD on `A` will not produce a singular value
1834 exactly equal to 0 in general: any difference of the smallest SVD value from
1835 0 will be caused by numerical imprecision in the calculation of the SVD.
1836 Our threshold for small SVD values takes this numerical imprecision into
1837 account, and the default threshold will detect such numerical rank
1838 deficiency. The threshold may declare a matrix `A` rank deficient even if
1839 the linear combination of some columns of `A` is not exactly equal to
1840 another column of `A` but only numerically very close to another column of
1841 `A`.
1842
1843 We chose our default threshold because it is in wide use. Other thresholds
1844 are possible. For example, elsewhere in the 2007 edition of *Numerical
1845 recipes* there is an alternative threshold of ``S.max() *
1846 np.finfo(A.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe
1847 this threshold as being based on "expected roundoff error" (p 71).
1848
1849 The thresholds above deal with floating point roundoff error in the
1850 calculation of the SVD. However, you may have more information about the
1851 sources of error in `A` that would make you consider other tolerance values
1852 to detect *effective* rank deficiency. The most useful measure of the
1853 tolerance depends on the operations you intend to use on your matrix. For
1854 example, if your data come from uncertain measurements with uncertainties
1855 greater than floating point epsilon, choosing a tolerance near that
1856 uncertainty may be preferable. The tolerance may be absolute if the
1857 uncertainties are absolute rather than relative.
1858
1859 References
1860 ----------
1861 .. [1] MATLAB reference documentation, "Rank"
1862 https://www.mathworks.com/help/techdoc/ref/rank.html
1863 .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
1864 "Numerical Recipes (3rd edition)", Cambridge University Press, 2007,
1865 page 795.
1866
1867 Examples
1868 --------
1869 >>> from numpy.linalg import matrix_rank
1870 >>> matrix_rank(np.eye(4)) # Full rank matrix
1871 4
1872 >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
1873 >>> matrix_rank(I)
1874 3
1875 >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
1876 1
1877 >>> matrix_rank(np.zeros((4,)))
1878 0
1879 """
1880 A = asarray(A)
1881 if A.ndim < 2:
1882 return int(not all(A==0))
1883 S = svd(A, compute_uv=False, hermitian=hermitian)
1884 if tol is None:
1885 tol = S.max(axis=-1, keepdims=True) * max(A.shape[-2:]) * finfo(S.dtype).eps
1886 else:
1887 tol = asarray(tol)[..., newaxis]
1888 return count_nonzero(S > tol, axis=-1)
1889
1890
1891# Generalized inverse
1892
1893def _pinv_dispatcher(a, rcond=None, hermitian=None):
1894 return (a,)
1895
1896
1897@array_function_dispatch(_pinv_dispatcher)
1898def pinv(a, rcond=1e-15, hermitian=False):
1899 """
1900 Compute the (Moore-Penrose) pseudo-inverse of a matrix.
1901
1902 Calculate the generalized inverse of a matrix using its
1903 singular-value decomposition (SVD) and including all
1904 *large* singular values.
1905
1906 .. versionchanged:: 1.14
1907 Can now operate on stacks of matrices
1908
1909 Parameters
1910 ----------
1911 a : (..., M, N) array_like
1912 Matrix or stack of matrices to be pseudo-inverted.
1913 rcond : (...) array_like of float
1914 Cutoff for small singular values.
1915 Singular values less than or equal to
1916 ``rcond * largest_singular_value`` are set to zero.
1917 Broadcasts against the stack of matrices.
1918 hermitian : bool, optional
1919 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1920 enabling a more efficient method for finding singular values.
1921 Defaults to False.
1922
1923 .. versionadded:: 1.17.0
1924
1925 Returns
1926 -------
1927 B : (..., N, M) ndarray
1928 The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
1929 is `B`.
1930
1931 Raises
1932 ------
1933 LinAlgError
1934 If the SVD computation does not converge.
1935
1936 See Also
1937 --------
1938 scipy.linalg.pinv : Similar function in SciPy.
1939 scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
1940 Hermitian matrix.
1941
1942 Notes
1943 -----
1944 The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
1945 defined as: "the matrix that 'solves' [the least-squares problem]
1946 :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
1947 :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
1948
1949 It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
1950 value decomposition of A, then
1951 :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
1952 orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
1953 of A's so-called singular values, (followed, typically, by
1954 zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
1955 consisting of the reciprocals of A's singular values
1956 (again, followed by zeros). [1]_
1957
1958 References
1959 ----------
1960 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1961 FL, Academic Press, Inc., 1980, pp. 139-142.
1962
1963 Examples
1964 --------
1965 The following example checks that ``a * a+ * a == a`` and
1966 ``a+ * a * a+ == a+``:
1967
1968 >>> a = np.random.randn(9, 6)
1969 >>> B = np.linalg.pinv(a)
1970 >>> np.allclose(a, np.dot(a, np.dot(B, a)))
1971 True
1972 >>> np.allclose(B, np.dot(B, np.dot(a, B)))
1973 True
1974
1975 """
1976 a, wrap = _makearray(a)
1977 rcond = asarray(rcond)
1978 if _is_empty_2d(a):
1979 m, n = a.shape[-2:]
1980 res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
1981 return wrap(res)
1982 a = a.conjugate()
1983 u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
1984
1985 # discard small singular values
1986 cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True)
1987 large = s > cutoff
1988 s = divide(1, s, where=large, out=s)
1989 s[~large] = 0
1990
1991 res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u)))
1992 return wrap(res)
1993
1994
1995# Determinant
1996
1997
1998@array_function_dispatch(_unary_dispatcher)
1999def slogdet(a):
2000 """
2001 Compute the sign and (natural) logarithm of the determinant of an array.
2002
2003 If an array has a very small or very large determinant, then a call to
2004 `det` may overflow or underflow. This routine is more robust against such
2005 issues, because it computes the logarithm of the determinant rather than
2006 the determinant itself.
2007
2008 Parameters
2009 ----------
2010 a : (..., M, M) array_like
2011 Input array, has to be a square 2-D array.
2012
2013 Returns
2014 -------
2015 sign : (...) array_like
2016 A number representing the sign of the determinant. For a real matrix,
2017 this is 1, 0, or -1. For a complex matrix, this is a complex number
2018 with absolute value 1 (i.e., it is on the unit circle), or else 0.
2019 logdet : (...) array_like
2020 The natural log of the absolute value of the determinant.
2021
2022 If the determinant is zero, then `sign` will be 0 and `logdet` will be
2023 -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``.
2024
2025 See Also
2026 --------
2027 det
2028
2029 Notes
2030 -----
2031
2032 .. versionadded:: 1.8.0
2033
2034 Broadcasting rules apply, see the `numpy.linalg` documentation for
2035 details.
2036
2037 .. versionadded:: 1.6.0
2038
2039 The determinant is computed via LU factorization using the LAPACK
2040 routine ``z/dgetrf``.
2041
2042
2043 Examples
2044 --------
2045 The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``:
2046
2047 >>> a = np.array([[1, 2], [3, 4]])
2048 >>> (sign, logdet) = np.linalg.slogdet(a)
2049 >>> (sign, logdet)
2050 (-1, 0.69314718055994529) # may vary
2051 >>> sign * np.exp(logdet)
2052 -2.0
2053
2054 Computing log-determinants for a stack of matrices:
2055
2056 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2057 >>> a.shape
2058 (3, 2, 2)
2059 >>> sign, logdet = np.linalg.slogdet(a)
2060 >>> (sign, logdet)
2061 (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154]))
2062 >>> sign * np.exp(logdet)
2063 array([-2., -3., -8.])
2064
2065 This routine succeeds where ordinary `det` does not:
2066
2067 >>> np.linalg.det(np.eye(500) * 0.1)
2068 0.0
2069 >>> np.linalg.slogdet(np.eye(500) * 0.1)
2070 (1, -1151.2925464970228)
2071
2072 """
2073 a = asarray(a)
2074 _assert_stacked_2d(a)
2075 _assert_stacked_square(a)
2076 t, result_t = _commonType(a)
2077 real_t = _realType(result_t)
2078 signature = 'D->Dd' if isComplexType(t) else 'd->dd'
2079 sign, logdet = _umath_linalg.slogdet(a, signature=signature)
2080 sign = sign.astype(result_t, copy=False)
2081 logdet = logdet.astype(real_t, copy=False)
2082 return sign, logdet
2083
2084
2085@array_function_dispatch(_unary_dispatcher)
2086def det(a):
2087 """
2088 Compute the determinant of an array.
2089
2090 Parameters
2091 ----------
2092 a : (..., M, M) array_like
2093 Input array to compute determinants for.
2094
2095 Returns
2096 -------
2097 det : (...) array_like
2098 Determinant of `a`.
2099
2100 See Also
2101 --------
2102 slogdet : Another way to represent the determinant, more suitable
2103 for large matrices where underflow/overflow may occur.
2104 scipy.linalg.det : Similar function in SciPy.
2105
2106 Notes
2107 -----
2108
2109 .. versionadded:: 1.8.0
2110
2111 Broadcasting rules apply, see the `numpy.linalg` documentation for
2112 details.
2113
2114 The determinant is computed via LU factorization using the LAPACK
2115 routine ``z/dgetrf``.
2116
2117 Examples
2118 --------
2119 The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
2120
2121 >>> a = np.array([[1, 2], [3, 4]])
2122 >>> np.linalg.det(a)
2123 -2.0 # may vary
2124
2125 Computing determinants for a stack of matrices:
2126
2127 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2128 >>> a.shape
2129 (3, 2, 2)
2130 >>> np.linalg.det(a)
2131 array([-2., -3., -8.])
2132
2133 """
2134 a = asarray(a)
2135 _assert_stacked_2d(a)
2136 _assert_stacked_square(a)
2137 t, result_t = _commonType(a)
2138 signature = 'D->D' if isComplexType(t) else 'd->d'
2139 r = _umath_linalg.det(a, signature=signature)
2140 r = r.astype(result_t, copy=False)
2141 return r
2142
2143
2144# Linear Least Squares
2145
2146def _lstsq_dispatcher(a, b, rcond=None):
2147 return (a, b)
2148
2149
2150@array_function_dispatch(_lstsq_dispatcher)
2151def lstsq(a, b, rcond="warn"):
2152 r"""
2153 Return the least-squares solution to a linear matrix equation.
2154
2155 Computes the vector `x` that approximately solves the equation
2156 ``a @ x = b``. The equation may be under-, well-, or over-determined
2157 (i.e., the number of linearly independent rows of `a` can be less than,
2158 equal to, or greater than its number of linearly independent columns).
2159 If `a` is square and of full rank, then `x` (but for round-off error)
2160 is the "exact" solution of the equation. Else, `x` minimizes the
2161 Euclidean 2-norm :math:`||b - ax||`. If there are multiple minimizing
2162 solutions, the one with the smallest 2-norm :math:`||x||` is returned.
2163
2164 Parameters
2165 ----------
2166 a : (M, N) array_like
2167 "Coefficient" matrix.
2168 b : {(M,), (M, K)} array_like
2169 Ordinate or "dependent variable" values. If `b` is two-dimensional,
2170 the least-squares solution is calculated for each of the `K` columns
2171 of `b`.
2172 rcond : float, optional
2173 Cut-off ratio for small singular values of `a`.
2174 For the purposes of rank determination, singular values are treated
2175 as zero if they are smaller than `rcond` times the largest singular
2176 value of `a`.
2177
2178 .. versionchanged:: 1.14.0
2179 If not set, a FutureWarning is given. The previous default
2180 of ``-1`` will use the machine precision as `rcond` parameter,
2181 the new default will use the machine precision times `max(M, N)`.
2182 To silence the warning and use the new default, use ``rcond=None``,
2183 to keep using the old behavior, use ``rcond=-1``.
2184
2185 Returns
2186 -------
2187 x : {(N,), (N, K)} ndarray
2188 Least-squares solution. If `b` is two-dimensional,
2189 the solutions are in the `K` columns of `x`.
2190 residuals : {(1,), (K,), (0,)} ndarray
2191 Sums of squared residuals: Squared Euclidean 2-norm for each column in
2192 ``b - a @ x``.
2193 If the rank of `a` is < N or M <= N, this is an empty array.
2194 If `b` is 1-dimensional, this is a (1,) shape array.
2195 Otherwise the shape is (K,).
2196 rank : int
2197 Rank of matrix `a`.
2198 s : (min(M, N),) ndarray
2199 Singular values of `a`.
2200
2201 Raises
2202 ------
2203 LinAlgError
2204 If computation does not converge.
2205
2206 See Also
2207 --------
2208 scipy.linalg.lstsq : Similar function in SciPy.
2209
2210 Notes
2211 -----
2212 If `b` is a matrix, then all array results are returned as matrices.
2213
2214 Examples
2215 --------
2216 Fit a line, ``y = mx + c``, through some noisy data-points:
2217
2218 >>> x = np.array([0, 1, 2, 3])
2219 >>> y = np.array([-1, 0.2, 0.9, 2.1])
2220
2221 By examining the coefficients, we see that the line should have a
2222 gradient of roughly 1 and cut the y-axis at, more or less, -1.
2223
2224 We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]``
2225 and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`:
2226
2227 >>> A = np.vstack([x, np.ones(len(x))]).T
2228 >>> A
2229 array([[ 0., 1.],
2230 [ 1., 1.],
2231 [ 2., 1.],
2232 [ 3., 1.]])
2233
2234 >>> m, c = np.linalg.lstsq(A, y, rcond=None)[0]
2235 >>> m, c
2236 (1.0 -0.95) # may vary
2237
2238 Plot the data along with the fitted line:
2239
2240 >>> import matplotlib.pyplot as plt
2241 >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10)
2242 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line')
2243 >>> _ = plt.legend()
2244 >>> plt.show()
2245
2246 """
2247 a, _ = _makearray(a)
2248 b, wrap = _makearray(b)
2249 is_1d = b.ndim == 1
2250 if is_1d:
2251 b = b[:, newaxis]
2252 _assert_2d(a, b)
2253 m, n = a.shape[-2:]
2254 m2, n_rhs = b.shape[-2:]
2255 if m != m2:
2256 raise LinAlgError('Incompatible dimensions')
2257
2258 t, result_t = _commonType(a, b)
2259 result_real_t = _realType(result_t)
2260
2261 # Determine default rcond value
2262 if rcond == "warn":
2263 # 2017-08-19, 1.14.0
2264 warnings.warn("`rcond` parameter will change to the default of "
2265 "machine precision times ``max(M, N)`` where M and N "
2266 "are the input matrix dimensions.\n"
2267 "To use the future default and silence this warning "
2268 "we advise to pass `rcond=None`, to keep using the old, "
2269 "explicitly pass `rcond=-1`.",
2270 FutureWarning, stacklevel=3)
2271 rcond = -1
2272 if rcond is None:
2273 rcond = finfo(t).eps * max(n, m)
2274
2275 if m <= n:
2276 gufunc = _umath_linalg.lstsq_m
2277 else:
2278 gufunc = _umath_linalg.lstsq_n
2279
2280 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid'
2281 extobj = get_linalg_error_extobj(_raise_linalgerror_lstsq)
2282 if n_rhs == 0:
2283 # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis
2284 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
2285 x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
2286 if m == 0:
2287 x[...] = 0
2288 if n_rhs == 0:
2289 # remove the item we added
2290 x = x[..., :n_rhs]
2291 resids = resids[..., :n_rhs]
2292
2293 # remove the axis we added
2294 if is_1d:
2295 x = x.squeeze(axis=-1)
2296 # we probably should squeeze resids too, but we can't
2297 # without breaking compatibility.
2298
2299 # as documented
2300 if rank != n or m <= n:
2301 resids = array([], result_real_t)
2302
2303 # coerce output arrays
2304 s = s.astype(result_real_t, copy=False)
2305 resids = resids.astype(result_real_t, copy=False)
2306 x = x.astype(result_t, copy=True) # Copying lets the memory in r_parts be freed
2307 return wrap(x), wrap(resids), rank, s
2308
2309
2310def _multi_svd_norm(x, row_axis, col_axis, op):
2311 """Compute a function of the singular values of the 2-D matrices in `x`.
2312
2313 This is a private utility function used by `numpy.linalg.norm()`.
2314
2315 Parameters
2316 ----------
2317 x : ndarray
2318 row_axis, col_axis : int
2319 The axes of `x` that hold the 2-D matrices.
2320 op : callable
2321 This should be either numpy.amin or `numpy.amax` or `numpy.sum`.
2322
2323 Returns
2324 -------
2325 result : float or ndarray
2326 If `x` is 2-D, the return values is a float.
2327 Otherwise, it is an array with ``x.ndim - 2`` dimensions.
2328 The return values are either the minimum or maximum or sum of the
2329 singular values of the matrices, depending on whether `op`
2330 is `numpy.amin` or `numpy.amax` or `numpy.sum`.
2331
2332 """
2333 y = moveaxis(x, (row_axis, col_axis), (-2, -1))
2334 result = op(svd(y, compute_uv=False), axis=-1)
2335 return result
2336
2337
2338def _norm_dispatcher(x, ord=None, axis=None, keepdims=None):
2339 return (x,)
2340
2341
2342@array_function_dispatch(_norm_dispatcher)
2343def norm(x, ord=None, axis=None, keepdims=False):
2344 """
2345 Matrix or vector norm.
2346
2347 This function is able to return one of eight different matrix norms,
2348 or one of an infinite number of vector norms (described below), depending
2349 on the value of the ``ord`` parameter.
2350
2351 Parameters
2352 ----------
2353 x : array_like
2354 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
2355 is None. If both `axis` and `ord` are None, the 2-norm of
2356 ``x.ravel`` will be returned.
2357 ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
2358 Order of the norm (see table under ``Notes``). inf means numpy's
2359 `inf` object. The default is None.
2360 axis : {None, int, 2-tuple of ints}, optional.
2361 If `axis` is an integer, it specifies the axis of `x` along which to
2362 compute the vector norms. If `axis` is a 2-tuple, it specifies the
2363 axes that hold 2-D matrices, and the matrix norms of these matrices
2364 are computed. If `axis` is None then either a vector norm (when `x`
2365 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
2366 is None.
2367
2368 .. versionadded:: 1.8.0
2369
2370 keepdims : bool, optional
2371 If this is set to True, the axes which are normed over are left in the
2372 result as dimensions with size one. With this option the result will
2373 broadcast correctly against the original `x`.
2374
2375 .. versionadded:: 1.10.0
2376
2377 Returns
2378 -------
2379 n : float or ndarray
2380 Norm of the matrix or vector(s).
2381
2382 See Also
2383 --------
2384 scipy.linalg.norm : Similar function in SciPy.
2385
2386 Notes
2387 -----
2388 For values of ``ord < 1``, the result is, strictly speaking, not a
2389 mathematical 'norm', but it may still be useful for various numerical
2390 purposes.
2391
2392 The following norms can be calculated:
2393
2394 ===== ============================ ==========================
2395 ord norm for matrices norm for vectors
2396 ===== ============================ ==========================
2397 None Frobenius norm 2-norm
2398 'fro' Frobenius norm --
2399 'nuc' nuclear norm --
2400 inf max(sum(abs(x), axis=1)) max(abs(x))
2401 -inf min(sum(abs(x), axis=1)) min(abs(x))
2402 0 -- sum(x != 0)
2403 1 max(sum(abs(x), axis=0)) as below
2404 -1 min(sum(abs(x), axis=0)) as below
2405 2 2-norm (largest sing. value) as below
2406 -2 smallest singular value as below
2407 other -- sum(abs(x)**ord)**(1./ord)
2408 ===== ============================ ==========================
2409
2410 The Frobenius norm is given by [1]_:
2411
2412 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
2413
2414 The nuclear norm is the sum of the singular values.
2415
2416 Both the Frobenius and nuclear norm orders are only defined for
2417 matrices and raise a ValueError when ``x.ndim != 2``.
2418
2419 References
2420 ----------
2421 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
2422 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
2423
2424 Examples
2425 --------
2426 >>> from numpy import linalg as LA
2427 >>> a = np.arange(9) - 4
2428 >>> a
2429 array([-4, -3, -2, ..., 2, 3, 4])
2430 >>> b = a.reshape((3, 3))
2431 >>> b
2432 array([[-4, -3, -2],
2433 [-1, 0, 1],
2434 [ 2, 3, 4]])
2435
2436 >>> LA.norm(a)
2437 7.745966692414834
2438 >>> LA.norm(b)
2439 7.745966692414834
2440 >>> LA.norm(b, 'fro')
2441 7.745966692414834
2442 >>> LA.norm(a, np.inf)
2443 4.0
2444 >>> LA.norm(b, np.inf)
2445 9.0
2446 >>> LA.norm(a, -np.inf)
2447 0.0
2448 >>> LA.norm(b, -np.inf)
2449 2.0
2450
2451 >>> LA.norm(a, 1)
2452 20.0
2453 >>> LA.norm(b, 1)
2454 7.0
2455 >>> LA.norm(a, -1)
2456 -4.6566128774142013e-010
2457 >>> LA.norm(b, -1)
2458 6.0
2459 >>> LA.norm(a, 2)
2460 7.745966692414834
2461 >>> LA.norm(b, 2)
2462 7.3484692283495345
2463
2464 >>> LA.norm(a, -2)
2465 0.0
2466 >>> LA.norm(b, -2)
2467 1.8570331885190563e-016 # may vary
2468 >>> LA.norm(a, 3)
2469 5.8480354764257312 # may vary
2470 >>> LA.norm(a, -3)
2471 0.0
2472
2473 Using the `axis` argument to compute vector norms:
2474
2475 >>> c = np.array([[ 1, 2, 3],
2476 ... [-1, 1, 4]])
2477 >>> LA.norm(c, axis=0)
2478 array([ 1.41421356, 2.23606798, 5. ])
2479 >>> LA.norm(c, axis=1)
2480 array([ 3.74165739, 4.24264069])
2481 >>> LA.norm(c, ord=1, axis=1)
2482 array([ 6., 6.])
2483
2484 Using the `axis` argument to compute matrix norms:
2485
2486 >>> m = np.arange(8).reshape(2,2,2)
2487 >>> LA.norm(m, axis=(1,2))
2488 array([ 3.74165739, 11.22497216])
2489 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :])
2490 (3.7416573867739413, 11.224972160321824)
2491
2492 """
2493 x = asarray(x)
2494
2495 if not issubclass(x.dtype.type, (inexact, object_)):
2496 x = x.astype(float)
2497
2498 # Immediately handle some default, simple, fast, and common cases.
2499 if axis is None:
2500 ndim = x.ndim
2501 if ((ord is None) or
2502 (ord in ('f', 'fro') and ndim == 2) or
2503 (ord == 2 and ndim == 1)):
2504
2505 x = x.ravel(order='K')
2506 if isComplexType(x.dtype.type):
2507 x_real = x.real
2508 x_imag = x.imag
2509 sqnorm = x_real.dot(x_real) + x_imag.dot(x_imag)
2510 else:
2511 sqnorm = x.dot(x)
2512 ret = sqrt(sqnorm)
2513 if keepdims:
2514 ret = ret.reshape(ndim*[1])
2515 return ret
2516
2517 # Normalize the `axis` argument to a tuple.
2518 nd = x.ndim
2519 if axis is None:
2520 axis = tuple(range(nd))
2521 elif not isinstance(axis, tuple):
2522 try:
2523 axis = int(axis)
2524 except Exception as e:
2525 raise TypeError("'axis' must be None, an integer or a tuple of integers") from e
2526 axis = (axis,)
2527
2528 if len(axis) == 1:
2529 if ord == Inf:
2530 return abs(x).max(axis=axis, keepdims=keepdims)
2531 elif ord == -Inf:
2532 return abs(x).min(axis=axis, keepdims=keepdims)
2533 elif ord == 0:
2534 # Zero norm
2535 return (x != 0).astype(x.real.dtype).sum(axis=axis, keepdims=keepdims)
2536 elif ord == 1:
2537 # special case for speedup
2538 return add.reduce(abs(x), axis=axis, keepdims=keepdims)
2539 elif ord is None or ord == 2:
2540 # special case for speedup
2541 s = (x.conj() * x).real
2542 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
2543 # None of the str-type keywords for ord ('fro', 'nuc')
2544 # are valid for vectors
2545 elif isinstance(ord, str):
2546 raise ValueError(f"Invalid norm order '{ord}' for vectors")
2547 else:
2548 absx = abs(x)
2549 absx **= ord
2550 ret = add.reduce(absx, axis=axis, keepdims=keepdims)
2551 ret **= reciprocal(ord, dtype=ret.dtype)
2552 return ret
2553 elif len(axis) == 2:
2554 row_axis, col_axis = axis
2555 row_axis = normalize_axis_index(row_axis, nd)
2556 col_axis = normalize_axis_index(col_axis, nd)
2557 if row_axis == col_axis:
2558 raise ValueError('Duplicate axes given.')
2559 if ord == 2:
2560 ret = _multi_svd_norm(x, row_axis, col_axis, amax)
2561 elif ord == -2:
2562 ret = _multi_svd_norm(x, row_axis, col_axis, amin)
2563 elif ord == 1:
2564 if col_axis > row_axis:
2565 col_axis -= 1
2566 ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis)
2567 elif ord == Inf:
2568 if row_axis > col_axis:
2569 row_axis -= 1
2570 ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis)
2571 elif ord == -1:
2572 if col_axis > row_axis:
2573 col_axis -= 1
2574 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis)
2575 elif ord == -Inf:
2576 if row_axis > col_axis:
2577 row_axis -= 1
2578 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis)
2579 elif ord in [None, 'fro', 'f']:
2580 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis))
2581 elif ord == 'nuc':
2582 ret = _multi_svd_norm(x, row_axis, col_axis, sum)
2583 else:
2584 raise ValueError("Invalid norm order for matrices.")
2585 if keepdims:
2586 ret_shape = list(x.shape)
2587 ret_shape[axis[0]] = 1
2588 ret_shape[axis[1]] = 1
2589 ret = ret.reshape(ret_shape)
2590 return ret
2591 else:
2592 raise ValueError("Improper number of dimensions to norm.")
2593
2594
2595# multi_dot
2596
2597def _multidot_dispatcher(arrays, *, out=None):
2598 yield from arrays
2599 yield out
2600
2601
2602@array_function_dispatch(_multidot_dispatcher)
2603def multi_dot(arrays, *, out=None):
2604 """
2605 Compute the dot product of two or more arrays in a single function call,
2606 while automatically selecting the fastest evaluation order.
2607
2608 `multi_dot` chains `numpy.dot` and uses optimal parenthesization
2609 of the matrices [1]_ [2]_. Depending on the shapes of the matrices,
2610 this can speed up the multiplication a lot.
2611
2612 If the first argument is 1-D it is treated as a row vector.
2613 If the last argument is 1-D it is treated as a column vector.
2614 The other arguments must be 2-D.
2615
2616 Think of `multi_dot` as::
2617
2618 def multi_dot(arrays): return functools.reduce(np.dot, arrays)
2619
2620
2621 Parameters
2622 ----------
2623 arrays : sequence of array_like
2624 If the first argument is 1-D it is treated as row vector.
2625 If the last argument is 1-D it is treated as column vector.
2626 The other arguments must be 2-D.
2627 out : ndarray, optional
2628 Output argument. This must have the exact kind that would be returned
2629 if it was not used. In particular, it must have the right type, must be
2630 C-contiguous, and its dtype must be the dtype that would be returned
2631 for `dot(a, b)`. This is a performance feature. Therefore, if these
2632 conditions are not met, an exception is raised, instead of attempting
2633 to be flexible.
2634
2635 .. versionadded:: 1.19.0
2636
2637 Returns
2638 -------
2639 output : ndarray
2640 Returns the dot product of the supplied arrays.
2641
2642 See Also
2643 --------
2644 numpy.dot : dot multiplication with two arguments.
2645
2646 References
2647 ----------
2648
2649 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378
2650 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication
2651
2652 Examples
2653 --------
2654 `multi_dot` allows you to write::
2655
2656 >>> from numpy.linalg import multi_dot
2657 >>> # Prepare some data
2658 >>> A = np.random.random((10000, 100))
2659 >>> B = np.random.random((100, 1000))
2660 >>> C = np.random.random((1000, 5))
2661 >>> D = np.random.random((5, 333))
2662 >>> # the actual dot multiplication
2663 >>> _ = multi_dot([A, B, C, D])
2664
2665 instead of::
2666
2667 >>> _ = np.dot(np.dot(np.dot(A, B), C), D)
2668 >>> # or
2669 >>> _ = A.dot(B).dot(C).dot(D)
2670
2671 Notes
2672 -----
2673 The cost for a matrix multiplication can be calculated with the
2674 following function::
2675
2676 def cost(A, B):
2677 return A.shape[0] * A.shape[1] * B.shape[1]
2678
2679 Assume we have three matrices
2680 :math:`A_{10x100}, B_{100x5}, C_{5x50}`.
2681
2682 The costs for the two different parenthesizations are as follows::
2683
2684 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500
2685 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000
2686
2687 """
2688 n = len(arrays)
2689 # optimization only makes sense for len(arrays) > 2
2690 if n < 2:
2691 raise ValueError("Expecting at least two arrays.")
2692 elif n == 2:
2693 return dot(arrays[0], arrays[1], out=out)
2694
2695 arrays = [asanyarray(a) for a in arrays]
2696
2697 # save original ndim to reshape the result array into the proper form later
2698 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim
2699 # Explicitly convert vectors to 2D arrays to keep the logic of the internal
2700 # _multi_dot_* functions as simple as possible.
2701 if arrays[0].ndim == 1:
2702 arrays[0] = atleast_2d(arrays[0])
2703 if arrays[-1].ndim == 1:
2704 arrays[-1] = atleast_2d(arrays[-1]).T
2705 _assert_2d(*arrays)
2706
2707 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order
2708 if n == 3:
2709 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out)
2710 else:
2711 order = _multi_dot_matrix_chain_order(arrays)
2712 result = _multi_dot(arrays, order, 0, n - 1, out=out)
2713
2714 # return proper shape
2715 if ndim_first == 1 and ndim_last == 1:
2716 return result[0, 0] # scalar
2717 elif ndim_first == 1 or ndim_last == 1:
2718 return result.ravel() # 1-D
2719 else:
2720 return result
2721
2722
2723def _multi_dot_three(A, B, C, out=None):
2724 """
2725 Find the best order for three arrays and do the multiplication.
2726
2727 For three arguments `_multi_dot_three` is approximately 15 times faster
2728 than `_multi_dot_matrix_chain_order`
2729
2730 """
2731 a0, a1b0 = A.shape
2732 b1c0, c1 = C.shape
2733 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1
2734 cost1 = a0 * b1c0 * (a1b0 + c1)
2735 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1
2736 cost2 = a1b0 * c1 * (a0 + b1c0)
2737
2738 if cost1 < cost2:
2739 return dot(dot(A, B), C, out=out)
2740 else:
2741 return dot(A, dot(B, C), out=out)
2742
2743
2744def _multi_dot_matrix_chain_order(arrays, return_costs=False):
2745 """
2746 Return a np.array that encodes the optimal order of mutiplications.
2747
2748 The optimal order array is then used by `_multi_dot()` to do the
2749 multiplication.
2750
2751 Also return the cost matrix if `return_costs` is `True`
2752
2753 The implementation CLOSELY follows Cormen, "Introduction to Algorithms",
2754 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices.
2755
2756 cost[i, j] = min([
2757 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix)
2758 for k in range(i, j)])
2759
2760 """
2761 n = len(arrays)
2762 # p stores the dimensions of the matrices
2763 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50]
2764 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]]
2765 # m is a matrix of costs of the subproblems
2766 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j}
2767 m = zeros((n, n), dtype=double)
2768 # s is the actual ordering
2769 # s[i, j] is the value of k at which we split the product A_i..A_j
2770 s = empty((n, n), dtype=intp)
2771
2772 for l in range(1, n):
2773 for i in range(n - l):
2774 j = i + l
2775 m[i, j] = Inf
2776 for k in range(i, j):
2777 q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1]
2778 if q < m[i, j]:
2779 m[i, j] = q
2780 s[i, j] = k # Note that Cormen uses 1-based index
2781
2782 return (s, m) if return_costs else s
2783
2784
2785def _multi_dot(arrays, order, i, j, out=None):
2786 """Actually do the multiplication with the given order."""
2787 if i == j:
2788 # the initial call with non-None out should never get here
2789 assert out is None
2790
2791 return arrays[i]
2792 else:
2793 return dot(_multi_dot(arrays, order, i, order[i, j]),
2794 _multi_dot(arrays, order, order[i, j] + 1, j),
2795 out=out)