codekingpro/portable-devtools
114k
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', 'svdvals', 'eig', 'eigh', 'lstsq', 'norm', 'qr', 'cond',
15 'matrix_rank', 'LinAlgError', 'multi_dot', 'trace', 'diagonal',
16 'cross', 'outer', 'tensordot', 'matmul', 'matrix_transpose',
17 'matrix_norm', 'vector_norm', 'vecdot']
18
19import functools
20import operator
21import warnings
22from typing import Any, NamedTuple
23
24from numpy._core import (
25 abs,
26 add,
27 all,
28 amax,
29 amin,
30 argsort,
31 array,
32 asanyarray,
33 asarray,
34 atleast_2d,
35 cdouble,
36 complexfloating,
37 count_nonzero,
38 cross as _core_cross,
39 csingle,
40 diagonal as _core_diagonal,
41 divide,
42 dot,
43 double,
44 empty,
45 empty_like,
46 errstate,
47 finfo,
48 inexact,
49 inf,
50 intc,
51 intp,
52 isfinite,
53 isnan,
54 matmul as _core_matmul,
55 matrix_transpose as _core_matrix_transpose,
56 moveaxis,
57 multiply,
58 newaxis,
59 object_,
60 outer as _core_outer,
61 overrides,
62 prod,
63 reciprocal,
64 sign,
65 single,
66 sort,
67 sqrt,
68 sum,
69 swapaxes,
70 tensordot as _core_tensordot,
71 trace as _core_trace,
72 transpose as _core_transpose,
73 vecdot as _core_vecdot,
74 zeros,
75)
76from numpy._globals import _NoValue
77from numpy._typing import NDArray
78from numpy._utils import set_module
79from numpy.lib._twodim_base_impl import eye, triu
80from numpy.lib.array_utils import normalize_axis_index, normalize_axis_tuple
81from numpy.linalg import _umath_linalg
82
83
84class EigResult(NamedTuple):
85 eigenvalues: NDArray[Any]
86 eigenvectors: NDArray[Any]
87
88class EighResult(NamedTuple):
89 eigenvalues: NDArray[Any]
90 eigenvectors: NDArray[Any]
91
92class QRResult(NamedTuple):
93 Q: NDArray[Any]
94 R: NDArray[Any]
95
96class SlogdetResult(NamedTuple):
97 sign: NDArray[Any]
98 logabsdet: NDArray[Any]
99
100class SVDResult(NamedTuple):
101 U: NDArray[Any]
102 S: NDArray[Any]
103 Vh: NDArray[Any]
104
105
106array_function_dispatch = functools.partial(
107 overrides.array_function_dispatch, module='numpy.linalg'
108)
109
110
111fortran_int = intc
112
113
114@set_module('numpy.linalg')
115class LinAlgError(ValueError):
116 """
117 Generic Python-exception-derived object raised by linalg functions.
118
119 General purpose exception class, derived from Python's ValueError
120 class, programmatically raised in linalg functions when a Linear
121 Algebra-related condition would prevent further correct execution of the
122 function.
123
124 Parameters
125 ----------
126 None
127
128 Examples
129 --------
130 >>> from numpy import linalg as LA
131 >>> LA.inv(np.zeros((2,2)))
132 Traceback (most recent call last):
133 File "<stdin>", line 1, in <module>
134 File "...linalg.py", line 350,
135 in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
136 File "...linalg.py", line 249,
137 in solve
138 raise LinAlgError('Singular matrix')
139 numpy.linalg.LinAlgError: Singular matrix
140
141 """
142
143
144def _raise_linalgerror_singular(err, flag):
145 raise LinAlgError("Singular matrix")
146
147def _raise_linalgerror_nonposdef(err, flag):
148 raise LinAlgError("Matrix is not positive definite")
149
150def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
151 raise LinAlgError("Eigenvalues did not converge")
152
153def _raise_linalgerror_svd_nonconvergence(err, flag):
154 raise LinAlgError("SVD did not converge")
155
156def _raise_linalgerror_lstsq(err, flag):
157 raise LinAlgError("SVD did not converge in Linear Least Squares")
158
159def _raise_linalgerror_qr(err, flag):
160 raise LinAlgError("Incorrect argument found while performing "
161 "QR factorization")
162
163
164def _makearray(a):
165 new = asarray(a)
166 wrap = getattr(a, "__array_wrap__", new.__array_wrap__)
167 return new, wrap
168
169def isComplexType(t):
170 return issubclass(t, complexfloating)
171
172
173_real_types_map = {single: single,
174 double: double,
175 csingle: single,
176 cdouble: double}
177
178_complex_types_map = {single: csingle,
179 double: cdouble,
180 csingle: csingle,
181 cdouble: cdouble}
182
183def _realType(t, default=double):
184 return _real_types_map.get(t, default)
185
186def _complexType(t, default=cdouble):
187 return _complex_types_map.get(t, default)
188
189def _commonType(*arrays):
190 # in lite version, use higher precision (always double or cdouble)
191 result_type = single
192 is_complex = False
193 for a in arrays:
194 type_ = a.dtype.type
195 if issubclass(type_, inexact):
196 if isComplexType(type_):
197 is_complex = True
198 rt = _realType(type_, default=None)
199 if rt is double:
200 result_type = double
201 elif rt is None:
202 # unsupported inexact scalar
203 raise TypeError(f"array type {a.dtype.name} is unsupported in linalg")
204 else:
205 result_type = double
206 if is_complex:
207 result_type = _complex_types_map[result_type]
208 return cdouble, result_type
209 else:
210 return double, result_type
211
212
213def _to_native_byte_order(*arrays):
214 ret = []
215 for arr in arrays:
216 if arr.dtype.byteorder not in ('=', '|'):
217 ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('=')))
218 else:
219 ret.append(arr)
220 if len(ret) == 1:
221 return ret[0]
222 else:
223 return ret
224
225
226def _assert_2d(*arrays):
227 for a in arrays:
228 if a.ndim != 2:
229 raise LinAlgError('%d-dimensional array given. Array must be '
230 'two-dimensional' % a.ndim)
231
232def _assert_stacked_2d(*arrays):
233 for a in arrays:
234 if a.ndim < 2:
235 raise LinAlgError('%d-dimensional array given. Array must be '
236 'at least two-dimensional' % a.ndim)
237
238def _assert_stacked_square(*arrays):
239 for a in arrays:
240 try:
241 m, n = a.shape[-2:]
242 except ValueError:
243 raise LinAlgError('%d-dimensional array given. Array must be '
244 'at least two-dimensional' % a.ndim)
245 if m != n:
246 raise LinAlgError('Last 2 dimensions of the array must be square')
247
248def _assert_finite(*arrays):
249 for a in arrays:
250 if not isfinite(a).all():
251 raise LinAlgError("Array must not contain infs or NaNs")
252
253def _is_empty_2d(arr):
254 # check size first for efficiency
255 return arr.size == 0 and prod(arr.shape[-2:]) == 0
256
257
258def transpose(a):
259 """
260 Transpose each matrix in a stack of matrices.
261
262 Unlike np.transpose, this only swaps the last two axes, rather than all of
263 them
264
265 Parameters
266 ----------
267 a : (...,M,N) array_like
268
269 Returns
270 -------
271 aT : (...,N,M) ndarray
272 """
273 return swapaxes(a, -1, -2)
274
275# Linear equations
276
277def _tensorsolve_dispatcher(a, b, axes=None):
278 return (a, b)
279
280
281@array_function_dispatch(_tensorsolve_dispatcher)
282def tensorsolve(a, b, axes=None):
283 """
284 Solve the tensor equation ``a x = b`` for x.
285
286 It is assumed that all indices of `x` are summed over in the product,
287 together with the rightmost indices of `a`, as is done in, for example,
288 ``tensordot(a, x, axes=x.ndim)``.
289
290 Parameters
291 ----------
292 a : array_like
293 Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals
294 the shape of that sub-tensor of `a` consisting of the appropriate
295 number of its rightmost indices, and must be such that
296 ``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be
297 'square').
298 b : array_like
299 Right-hand tensor, which can be of any shape.
300 axes : tuple of ints, optional
301 Axes in `a` to reorder to the right, before inversion.
302 If None (default), no reordering is done.
303
304 Returns
305 -------
306 x : ndarray, shape Q
307
308 Raises
309 ------
310 LinAlgError
311 If `a` is singular or not 'square' (in the above sense).
312
313 See Also
314 --------
315 numpy.tensordot, tensorinv, numpy.einsum
316
317 Examples
318 --------
319 >>> import numpy as np
320 >>> a = np.eye(2*3*4).reshape((2*3, 4, 2, 3, 4))
321 >>> rng = np.random.default_rng()
322 >>> b = rng.normal(size=(2*3, 4))
323 >>> x = np.linalg.tensorsolve(a, b)
324 >>> x.shape
325 (2, 3, 4)
326 >>> np.allclose(np.tensordot(a, x, axes=3), b)
327 True
328
329 """
330 a, wrap = _makearray(a)
331 b = asarray(b)
332 an = a.ndim
333
334 if axes is not None:
335 allaxes = list(range(an))
336 for k in axes:
337 allaxes.remove(k)
338 allaxes.insert(an, k)
339 a = a.transpose(allaxes)
340
341 oldshape = a.shape[-(an - b.ndim):]
342 prod = 1
343 for k in oldshape:
344 prod *= k
345
346 if a.size != prod ** 2:
347 raise LinAlgError(
348 "Input arrays must satisfy the requirement \
349 prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])"
350 )
351
352 a = a.reshape(prod, prod)
353 b = b.ravel()
354 res = wrap(solve(a, b))
355 res.shape = oldshape
356 return res
357
358
359def _solve_dispatcher(a, b):
360 return (a, b)
361
362
363@array_function_dispatch(_solve_dispatcher)
364def solve(a, b):
365 """
366 Solve a linear matrix equation, or system of linear scalar equations.
367
368 Computes the "exact" solution, `x`, of the well-determined, i.e., full
369 rank, linear matrix equation `ax = b`.
370
371 Parameters
372 ----------
373 a : (..., M, M) array_like
374 Coefficient matrix.
375 b : {(M,), (..., M, K)}, array_like
376 Ordinate or "dependent variable" values.
377
378 Returns
379 -------
380 x : {(..., M,), (..., M, K)} ndarray
381 Solution to the system a x = b. Returned shape is (..., M) if b is
382 shape (M,) and (..., M, K) if b is (..., M, K), where the "..." part is
383 broadcasted between a and b.
384
385 Raises
386 ------
387 LinAlgError
388 If `a` is singular or not square.
389
390 See Also
391 --------
392 scipy.linalg.solve : Similar function in SciPy.
393
394 Notes
395 -----
396 Broadcasting rules apply, see the `numpy.linalg` documentation for
397 details.
398
399 The solutions are computed using LAPACK routine ``_gesv``.
400
401 `a` must be square and of full-rank, i.e., all rows (or, equivalently,
402 columns) must be linearly independent; if either is not true, use
403 `lstsq` for the least-squares best "solution" of the
404 system/equation.
405
406 .. versionchanged:: 2.0
407
408 The b array is only treated as a shape (M,) column vector if it is
409 exactly 1-dimensional. In all other instances it is treated as a stack
410 of (M, K) matrices. Previously b would be treated as a stack of (M,)
411 vectors if b.ndim was equal to a.ndim - 1.
412
413 References
414 ----------
415 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
416 FL, Academic Press, Inc., 1980, pg. 22.
417
418 Examples
419 --------
420 Solve the system of equations:
421 ``x0 + 2 * x1 = 1`` and
422 ``3 * x0 + 5 * x1 = 2``:
423
424 >>> import numpy as np
425 >>> a = np.array([[1, 2], [3, 5]])
426 >>> b = np.array([1, 2])
427 >>> x = np.linalg.solve(a, b)
428 >>> x
429 array([-1., 1.])
430
431 Check that the solution is correct:
432
433 >>> np.allclose(np.dot(a, x), b)
434 True
435
436 """
437 a, _ = _makearray(a)
438 _assert_stacked_square(a)
439 b, wrap = _makearray(b)
440 t, result_t = _commonType(a, b)
441
442 # We use the b = (..., M,) logic, only if the number of extra dimensions
443 # match exactly
444 if b.ndim == 1:
445 gufunc = _umath_linalg.solve1
446 else:
447 gufunc = _umath_linalg.solve
448
449 signature = 'DD->D' if isComplexType(t) else 'dd->d'
450 with errstate(call=_raise_linalgerror_singular, invalid='call',
451 over='ignore', divide='ignore', under='ignore'):
452 r = gufunc(a, b, signature=signature)
453
454 return wrap(r.astype(result_t, copy=False))
455
456
457def _tensorinv_dispatcher(a, ind=None):
458 return (a,)
459
460
461@array_function_dispatch(_tensorinv_dispatcher)
462def tensorinv(a, ind=2):
463 """
464 Compute the 'inverse' of an N-dimensional array.
465
466 The result is an inverse for `a` relative to the tensordot operation
467 ``tensordot(a, b, ind)``, i. e., up to floating-point accuracy,
468 ``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the
469 tensordot operation.
470
471 Parameters
472 ----------
473 a : array_like
474 Tensor to 'invert'. Its shape must be 'square', i. e.,
475 ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.
476 ind : int, optional
477 Number of first indices that are involved in the inverse sum.
478 Must be a positive integer, default is 2.
479
480 Returns
481 -------
482 b : ndarray
483 `a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``.
484
485 Raises
486 ------
487 LinAlgError
488 If `a` is singular or not 'square' (in the above sense).
489
490 See Also
491 --------
492 numpy.tensordot, tensorsolve
493
494 Examples
495 --------
496 >>> import numpy as np
497 >>> a = np.eye(4*6).reshape((4, 6, 8, 3))
498 >>> ainv = np.linalg.tensorinv(a, ind=2)
499 >>> ainv.shape
500 (8, 3, 4, 6)
501 >>> rng = np.random.default_rng()
502 >>> b = rng.normal(size=(4, 6))
503 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
504 True
505
506 >>> a = np.eye(4*6).reshape((24, 8, 3))
507 >>> ainv = np.linalg.tensorinv(a, ind=1)
508 >>> ainv.shape
509 (8, 3, 24)
510 >>> rng = np.random.default_rng()
511 >>> b = rng.normal(size=24)
512 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
513 True
514
515 """
516 a = asarray(a)
517 oldshape = a.shape
518 prod = 1
519 if ind > 0:
520 invshape = oldshape[ind:] + oldshape[:ind]
521 for k in oldshape[ind:]:
522 prod *= k
523 else:
524 raise ValueError("Invalid ind argument.")
525 a = a.reshape(prod, -1)
526 ia = inv(a)
527 return ia.reshape(*invshape)
528
529
530# Matrix inversion
531
532def _unary_dispatcher(a):
533 return (a,)
534
535
536@array_function_dispatch(_unary_dispatcher)
537def inv(a):
538 """
539 Compute the inverse of a matrix.
540
541 Given a square matrix `a`, return the matrix `ainv` satisfying
542 ``a @ ainv = ainv @ a = eye(a.shape[0])``.
543
544 Parameters
545 ----------
546 a : (..., M, M) array_like
547 Matrix to be inverted.
548
549 Returns
550 -------
551 ainv : (..., M, M) ndarray or matrix
552 Inverse of the matrix `a`.
553
554 Raises
555 ------
556 LinAlgError
557 If `a` is not square or inversion fails.
558
559 See Also
560 --------
561 scipy.linalg.inv : Similar function in SciPy.
562 numpy.linalg.cond : Compute the condition number of a matrix.
563 numpy.linalg.svd : Compute the singular value decomposition of a matrix.
564
565 Notes
566 -----
567 Broadcasting rules apply, see the `numpy.linalg` documentation for
568 details.
569
570 If `a` is detected to be singular, a `LinAlgError` is raised. If `a` is
571 ill-conditioned, a `LinAlgError` may or may not be raised, and results may
572 be inaccurate due to floating-point errors.
573
574 References
575 ----------
576 .. [1] Wikipedia, "Condition number",
577 https://en.wikipedia.org/wiki/Condition_number
578
579 Examples
580 --------
581 >>> import numpy as np
582 >>> from numpy.linalg import inv
583 >>> a = np.array([[1., 2.], [3., 4.]])
584 >>> ainv = inv(a)
585 >>> np.allclose(a @ ainv, np.eye(2))
586 True
587 >>> np.allclose(ainv @ a, np.eye(2))
588 True
589
590 If a is a matrix object, then the return value is a matrix as well:
591
592 >>> ainv = inv(np.matrix(a))
593 >>> ainv
594 matrix([[-2. , 1. ],
595 [ 1.5, -0.5]])
596
597 Inverses of several matrices can be computed at once:
598
599 >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
600 >>> inv(a)
601 array([[[-2. , 1. ],
602 [ 1.5 , -0.5 ]],
603 [[-1.25, 0.75],
604 [ 0.75, -0.25]]])
605
606 If a matrix is close to singular, the computed inverse may not satisfy
607 ``a @ ainv = ainv @ a = eye(a.shape[0])`` even if a `LinAlgError`
608 is not raised:
609
610 >>> a = np.array([[2,4,6],[2,0,2],[6,8,14]])
611 >>> inv(a) # No errors raised
612 array([[-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],
613 [-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],
614 [ 1.12589991e+15, 5.62949953e+14, -5.62949953e+14]])
615 >>> a @ inv(a)
616 array([[ 0. , -0.5 , 0. ], # may vary
617 [-0.5 , 0.625, 0.25 ],
618 [ 0. , 0. , 1. ]])
619
620 To detect ill-conditioned matrices, you can use `numpy.linalg.cond` to
621 compute its *condition number* [1]_. The larger the condition number, the
622 more ill-conditioned the matrix is. As a rule of thumb, if the condition
623 number ``cond(a) = 10**k``, then you may lose up to ``k`` digits of
624 accuracy on top of what would be lost to the numerical method due to loss
625 of precision from arithmetic methods.
626
627 >>> from numpy.linalg import cond
628 >>> cond(a)
629 np.float64(8.659885634118668e+17) # may vary
630
631 It is also possible to detect ill-conditioning by inspecting the matrix's
632 singular values directly. The ratio between the largest and the smallest
633 singular value is the condition number:
634
635 >>> from numpy.linalg import svd
636 >>> sigma = svd(a, compute_uv=False) # Do not compute singular vectors
637 >>> sigma.max()/sigma.min()
638 8.659885634118668e+17 # may vary
639
640 """
641 a, wrap = _makearray(a)
642 _assert_stacked_square(a)
643 t, result_t = _commonType(a)
644
645 signature = 'D->D' if isComplexType(t) else 'd->d'
646 with errstate(call=_raise_linalgerror_singular, invalid='call',
647 over='ignore', divide='ignore', under='ignore'):
648 ainv = _umath_linalg.inv(a, signature=signature)
649 return wrap(ainv.astype(result_t, copy=False))
650
651
652def _matrix_power_dispatcher(a, n):
653 return (a,)
654
655
656@array_function_dispatch(_matrix_power_dispatcher)
657def matrix_power(a, n):
658 """
659 Raise a square matrix to the (integer) power `n`.
660
661 For positive integers `n`, the power is computed by repeated matrix
662 squarings and matrix multiplications. If ``n == 0``, the identity matrix
663 of the same shape as M is returned. If ``n < 0``, the inverse
664 is computed and then raised to the ``abs(n)``.
665
666 .. note:: Stacks of object matrices are not currently supported.
667
668 Parameters
669 ----------
670 a : (..., M, M) array_like
671 Matrix to be "powered".
672 n : int
673 The exponent can be any integer or long integer, positive,
674 negative, or zero.
675
676 Returns
677 -------
678 a**n : (..., M, M) ndarray or matrix object
679 The return value is the same shape and type as `M`;
680 if the exponent is positive or zero then the type of the
681 elements is the same as those of `M`. If the exponent is
682 negative the elements are floating-point.
683
684 Raises
685 ------
686 LinAlgError
687 For matrices that are not square or that (for negative powers) cannot
688 be inverted numerically.
689
690 Examples
691 --------
692 >>> import numpy as np
693 >>> from numpy.linalg import matrix_power
694 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
695 >>> matrix_power(i, 3) # should = -i
696 array([[ 0, -1],
697 [ 1, 0]])
698 >>> matrix_power(i, 0)
699 array([[1, 0],
700 [0, 1]])
701 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
702 array([[ 0., 1.],
703 [-1., 0.]])
704
705 Somewhat more sophisticated example
706
707 >>> q = np.zeros((4, 4))
708 >>> q[0:2, 0:2] = -i
709 >>> q[2:4, 2:4] = i
710 >>> q # one of the three quaternion units not equal to 1
711 array([[ 0., -1., 0., 0.],
712 [ 1., 0., 0., 0.],
713 [ 0., 0., 0., 1.],
714 [ 0., 0., -1., 0.]])
715 >>> matrix_power(q, 2) # = -np.eye(4)
716 array([[-1., 0., 0., 0.],
717 [ 0., -1., 0., 0.],
718 [ 0., 0., -1., 0.],
719 [ 0., 0., 0., -1.]])
720
721 """
722 a = asanyarray(a)
723 _assert_stacked_square(a)
724
725 try:
726 n = operator.index(n)
727 except TypeError as e:
728 raise TypeError("exponent must be an integer") from e
729
730 # Fall back on dot for object arrays. Object arrays are not supported by
731 # the current implementation of matmul using einsum
732 if a.dtype != object:
733 fmatmul = matmul
734 elif a.ndim == 2:
735 fmatmul = dot
736 else:
737 raise NotImplementedError(
738 "matrix_power not supported for stacks of object arrays")
739
740 if n == 0:
741 a = empty_like(a)
742 a[...] = eye(a.shape[-2], dtype=a.dtype)
743 return a
744
745 elif n < 0:
746 a = inv(a)
747 n = abs(n)
748
749 # short-cuts.
750 if n == 1:
751 return a
752
753 elif n == 2:
754 return fmatmul(a, a)
755
756 elif n == 3:
757 return fmatmul(fmatmul(a, a), a)
758
759 # Use binary decomposition to reduce the number of matrix multiplications.
760 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
761 # increasing powers of 2, and multiply into the result as needed.
762 z = result = None
763 while n > 0:
764 z = a if z is None else fmatmul(z, z)
765 n, bit = divmod(n, 2)
766 if bit:
767 result = z if result is None else fmatmul(result, z)
768
769 return result
770
771
772# Cholesky decomposition
773
774def _cholesky_dispatcher(a, /, *, upper=None):
775 return (a,)
776
777
778@array_function_dispatch(_cholesky_dispatcher)
779def cholesky(a, /, *, upper=False):
780 """
781 Cholesky decomposition.
782
783 Return the lower or upper Cholesky decomposition, ``L * L.H`` or
784 ``U.H * U``, of the square matrix ``a``, where ``L`` is lower-triangular,
785 ``U`` is upper-triangular, and ``.H`` is the conjugate transpose operator
786 (which is the ordinary transpose if ``a`` is real-valued). ``a`` must be
787 Hermitian (symmetric if real-valued) and positive-definite. No checking is
788 performed to verify whether ``a`` is Hermitian or not. In addition, only
789 the lower or upper-triangular and diagonal elements of ``a`` are used.
790 Only ``L`` or ``U`` is actually returned.
791
792 Parameters
793 ----------
794 a : (..., M, M) array_like
795 Hermitian (symmetric if all elements are real), positive-definite
796 input matrix.
797 upper : bool
798 If ``True``, the result must be the upper-triangular Cholesky factor.
799 If ``False``, the result must be the lower-triangular Cholesky factor.
800 Default: ``False``.
801
802 Returns
803 -------
804 L : (..., M, M) array_like
805 Lower or upper-triangular Cholesky factor of `a`. Returns a matrix
806 object if `a` is a matrix object.
807
808 Raises
809 ------
810 LinAlgError
811 If the decomposition fails, for example, if `a` is not
812 positive-definite.
813
814 See Also
815 --------
816 scipy.linalg.cholesky : Similar function in SciPy.
817 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
818 positive-definite matrix.
819 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
820 `scipy.linalg.cho_solve`.
821
822 Notes
823 -----
824 Broadcasting rules apply, see the `numpy.linalg` documentation for
825 details.
826
827 The Cholesky decomposition is often used as a fast way of solving
828
829 .. math:: A \\mathbf{x} = \\mathbf{b}
830
831 (when `A` is both Hermitian/symmetric and positive-definite).
832
833 First, we solve for :math:`\\mathbf{y}` in
834
835 .. math:: L \\mathbf{y} = \\mathbf{b},
836
837 and then for :math:`\\mathbf{x}` in
838
839 .. math:: L^{H} \\mathbf{x} = \\mathbf{y}.
840
841 Examples
842 --------
843 >>> import numpy as np
844 >>> A = np.array([[1,-2j],[2j,5]])
845 >>> A
846 array([[ 1.+0.j, -0.-2.j],
847 [ 0.+2.j, 5.+0.j]])
848 >>> L = np.linalg.cholesky(A)
849 >>> L
850 array([[1.+0.j, 0.+0.j],
851 [0.+2.j, 1.+0.j]])
852 >>> np.dot(L, L.T.conj()) # verify that L * L.H = A
853 array([[1.+0.j, 0.-2.j],
854 [0.+2.j, 5.+0.j]])
855 >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
856 >>> np.linalg.cholesky(A) # an ndarray object is returned
857 array([[1.+0.j, 0.+0.j],
858 [0.+2.j, 1.+0.j]])
859 >>> # But a matrix object is returned if A is a matrix object
860 >>> np.linalg.cholesky(np.matrix(A))
861 matrix([[ 1.+0.j, 0.+0.j],
862 [ 0.+2.j, 1.+0.j]])
863 >>> # The upper-triangular Cholesky factor can also be obtained.
864 >>> np.linalg.cholesky(A, upper=True)
865 array([[1.-0.j, 0.-2.j],
866 [0.-0.j, 1.-0.j]])
867
868 """
869 gufunc = _umath_linalg.cholesky_up if upper else _umath_linalg.cholesky_lo
870 a, wrap = _makearray(a)
871 _assert_stacked_square(a)
872 t, result_t = _commonType(a)
873 signature = 'D->D' if isComplexType(t) else 'd->d'
874 with errstate(call=_raise_linalgerror_nonposdef, invalid='call',
875 over='ignore', divide='ignore', under='ignore'):
876 r = gufunc(a, signature=signature)
877 return wrap(r.astype(result_t, copy=False))
878
879
880# outer product
881
882
883def _outer_dispatcher(x1, x2):
884 return (x1, x2)
885
886
887@array_function_dispatch(_outer_dispatcher)
888def outer(x1, x2, /):
889 """
890 Compute the outer product of two vectors.
891
892 This function is Array API compatible. Compared to ``np.outer``
893 it accepts 1-dimensional inputs only.
894
895 Parameters
896 ----------
897 x1 : (M,) array_like
898 One-dimensional input array of size ``N``.
899 Must have a numeric data type.
900 x2 : (N,) array_like
901 One-dimensional input array of size ``M``.
902 Must have a numeric data type.
903
904 Returns
905 -------
906 out : (M, N) ndarray
907 ``out[i, j] = a[i] * b[j]``
908
909 See also
910 --------
911 outer
912
913 Examples
914 --------
915 Make a (*very* coarse) grid for computing a Mandelbrot set:
916
917 >>> rl = np.linalg.outer(np.ones((5,)), np.linspace(-2, 2, 5))
918 >>> rl
919 array([[-2., -1., 0., 1., 2.],
920 [-2., -1., 0., 1., 2.],
921 [-2., -1., 0., 1., 2.],
922 [-2., -1., 0., 1., 2.],
923 [-2., -1., 0., 1., 2.]])
924 >>> im = np.linalg.outer(1j*np.linspace(2, -2, 5), np.ones((5,)))
925 >>> im
926 array([[0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j],
927 [0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j],
928 [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
929 [0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j],
930 [0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j]])
931 >>> grid = rl + im
932 >>> grid
933 array([[-2.+2.j, -1.+2.j, 0.+2.j, 1.+2.j, 2.+2.j],
934 [-2.+1.j, -1.+1.j, 0.+1.j, 1.+1.j, 2.+1.j],
935 [-2.+0.j, -1.+0.j, 0.+0.j, 1.+0.j, 2.+0.j],
936 [-2.-1.j, -1.-1.j, 0.-1.j, 1.-1.j, 2.-1.j],
937 [-2.-2.j, -1.-2.j, 0.-2.j, 1.-2.j, 2.-2.j]])
938
939 An example using a "vector" of letters:
940
941 >>> x = np.array(['a', 'b', 'c'], dtype=object)
942 >>> np.linalg.outer(x, [1, 2, 3])
943 array([['a', 'aa', 'aaa'],
944 ['b', 'bb', 'bbb'],
945 ['c', 'cc', 'ccc']], dtype=object)
946
947 """
948 x1 = asanyarray(x1)
949 x2 = asanyarray(x2)
950 if x1.ndim != 1 or x2.ndim != 1:
951 raise ValueError(
952 "Input arrays must be one-dimensional, but they are "
953 f"{x1.ndim=} and {x2.ndim=}."
954 )
955 return _core_outer(x1, x2, out=None)
956
957
958# QR decomposition
959
960
961def _qr_dispatcher(a, mode=None):
962 return (a,)
963
964
965@array_function_dispatch(_qr_dispatcher)
966def qr(a, mode='reduced'):
967 """
968 Compute the qr factorization of a matrix.
969
970 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is
971 upper-triangular.
972
973 Parameters
974 ----------
975 a : array_like, shape (..., M, N)
976 An array-like object with the dimensionality of at least 2.
977 mode : {'reduced', 'complete', 'r', 'raw'}, optional, default: 'reduced'
978 If K = min(M, N), then
979
980 * 'reduced' : returns Q, R with dimensions (..., M, K), (..., K, N)
981 * 'complete' : returns Q, R with dimensions (..., M, M), (..., M, N)
982 * 'r' : returns R only with dimensions (..., K, N)
983 * 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)
984
985 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
986 see the notes for more information. The default is 'reduced', and to
987 maintain backward compatibility with earlier versions of numpy both
988 it and the old default 'full' can be omitted. Note that array h
989 returned in 'raw' mode is transposed for calling Fortran. The
990 'economic' mode is deprecated. The modes 'full' and 'economic' may
991 be passed using only the first letter for backwards compatibility,
992 but all others must be spelled out. See the Notes for more
993 explanation.
994
995
996 Returns
997 -------
998 Q : ndarray of float or complex, optional
999 A matrix with orthonormal columns. When mode = 'complete' the
1000 result is an orthogonal/unitary matrix depending on whether or not
1001 a is real/complex. The determinant may be either +/- 1 in that
1002 case. In case the number of dimensions in the input array is
1003 greater than 2 then a stack of the matrices with above properties
1004 is returned.
1005 R : ndarray of float or complex, optional
1006 The upper-triangular matrix or a stack of upper-triangular
1007 matrices if the number of dimensions in the input array is greater
1008 than 2.
1009 (h, tau) : ndarrays of np.double or np.cdouble, optional
1010 The array h contains the Householder reflectors that generate q
1011 along with r. The tau array contains scaling factors for the
1012 reflectors. In the deprecated 'economic' mode only h is returned.
1013
1014 Raises
1015 ------
1016 LinAlgError
1017 If factoring fails.
1018
1019 See Also
1020 --------
1021 scipy.linalg.qr : Similar function in SciPy.
1022 scipy.linalg.rq : Compute RQ decomposition of a matrix.
1023
1024 Notes
1025 -----
1026 When mode is 'reduced' or 'complete', the result will be a namedtuple with
1027 the attributes ``Q`` and ``R``.
1028
1029 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
1030 ``dorgqr``, and ``zungqr``.
1031
1032 For more information on the qr factorization, see for example:
1033 https://en.wikipedia.org/wiki/QR_factorization
1034
1035 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
1036 `a` is of type `matrix`, all the return values will be matrices too.
1037
1038 New 'reduced', 'complete', and 'raw' options for mode were added in
1039 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In
1040 addition the options 'full' and 'economic' were deprecated. Because
1041 'full' was the previous default and 'reduced' is the new default,
1042 backward compatibility can be maintained by letting `mode` default.
1043 The 'raw' option was added so that LAPACK routines that can multiply
1044 arrays by q using the Householder reflectors can be used. Note that in
1045 this case the returned arrays are of type np.double or np.cdouble and
1046 the h array is transposed to be FORTRAN compatible. No routines using
1047 the 'raw' return are currently exposed by numpy, but some are available
1048 in lapack_lite and just await the necessary work.
1049
1050 Examples
1051 --------
1052 >>> import numpy as np
1053 >>> rng = np.random.default_rng()
1054 >>> a = rng.normal(size=(9, 6))
1055 >>> Q, R = np.linalg.qr(a)
1056 >>> np.allclose(a, np.dot(Q, R)) # a does equal QR
1057 True
1058 >>> R2 = np.linalg.qr(a, mode='r')
1059 >>> np.allclose(R, R2) # mode='r' returns the same R as mode='full'
1060 True
1061 >>> a = np.random.normal(size=(3, 2, 2)) # Stack of 2 x 2 matrices as input
1062 >>> Q, R = np.linalg.qr(a)
1063 >>> Q.shape
1064 (3, 2, 2)
1065 >>> R.shape
1066 (3, 2, 2)
1067 >>> np.allclose(a, np.matmul(Q, R))
1068 True
1069
1070 Example illustrating a common use of `qr`: solving of least squares
1071 problems
1072
1073 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for
1074 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
1075 and you'll see that it should be y0 = 0, m = 1.) The answer is provided
1076 by solving the over-determined matrix equation ``Ax = b``, where::
1077
1078 A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
1079 x = array([[y0], [m]])
1080 b = array([[1], [0], [2], [1]])
1081
1082 If A = QR such that Q is orthonormal (which is always possible via
1083 Gram-Schmidt), then ``x = inv(R) * (Q.T) * b``. (In numpy practice,
1084 however, we simply use `lstsq`.)
1085
1086 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
1087 >>> A
1088 array([[0, 1],
1089 [1, 1],
1090 [1, 1],
1091 [2, 1]])
1092 >>> b = np.array([1, 2, 2, 3])
1093 >>> Q, R = np.linalg.qr(A)
1094 >>> p = np.dot(Q.T, b)
1095 >>> np.dot(np.linalg.inv(R), p)
1096 array([ 1., 1.])
1097
1098 """
1099 if mode not in ('reduced', 'complete', 'r', 'raw'):
1100 if mode in ('f', 'full'):
1101 # 2013-04-01, 1.8
1102 msg = (
1103 "The 'full' option is deprecated in favor of 'reduced'.\n"
1104 "For backward compatibility let mode default."
1105 )
1106 warnings.warn(msg, DeprecationWarning, stacklevel=2)
1107 mode = 'reduced'
1108 elif mode in ('e', 'economic'):
1109 # 2013-04-01, 1.8
1110 msg = "The 'economic' option is deprecated."
1111 warnings.warn(msg, DeprecationWarning, stacklevel=2)
1112 mode = 'economic'
1113 else:
1114 raise ValueError(f"Unrecognized mode '{mode}'")
1115
1116 a, wrap = _makearray(a)
1117 _assert_stacked_2d(a)
1118 m, n = a.shape[-2:]
1119 t, result_t = _commonType(a)
1120 a = a.astype(t, copy=True)
1121 a = _to_native_byte_order(a)
1122 mn = min(m, n)
1123
1124 signature = 'D->D' if isComplexType(t) else 'd->d'
1125 with errstate(call=_raise_linalgerror_qr, invalid='call',
1126 over='ignore', divide='ignore', under='ignore'):
1127 tau = _umath_linalg.qr_r_raw(a, signature=signature)
1128
1129 # handle modes that don't return q
1130 if mode == 'r':
1131 r = triu(a[..., :mn, :])
1132 r = r.astype(result_t, copy=False)
1133 return wrap(r)
1134
1135 if mode == 'raw':
1136 q = transpose(a)
1137 q = q.astype(result_t, copy=False)
1138 tau = tau.astype(result_t, copy=False)
1139 return wrap(q), tau
1140
1141 if mode == 'economic':
1142 a = a.astype(result_t, copy=False)
1143 return wrap(a)
1144
1145 # mc is the number of columns in the resulting q
1146 # matrix. If the mode is complete then it is
1147 # same as number of rows, and if the mode is reduced,
1148 # then it is the minimum of number of rows and columns.
1149 if mode == 'complete' and m > n:
1150 mc = m
1151 gufunc = _umath_linalg.qr_complete
1152 else:
1153 mc = mn
1154 gufunc = _umath_linalg.qr_reduced
1155
1156 signature = 'DD->D' if isComplexType(t) else 'dd->d'
1157 with errstate(call=_raise_linalgerror_qr, invalid='call',
1158 over='ignore', divide='ignore', under='ignore'):
1159 q = gufunc(a, tau, signature=signature)
1160 r = triu(a[..., :mc, :])
1161
1162 q = q.astype(result_t, copy=False)
1163 r = r.astype(result_t, copy=False)
1164
1165 return QRResult(wrap(q), wrap(r))
1166
1167# Eigenvalues
1168
1169
1170@array_function_dispatch(_unary_dispatcher)
1171def eigvals(a):
1172 """
1173 Compute the eigenvalues of a general matrix.
1174
1175 Main difference between `eigvals` and `eig`: the eigenvectors aren't
1176 returned.
1177
1178 Parameters
1179 ----------
1180 a : (..., M, M) array_like
1181 A complex- or real-valued matrix whose eigenvalues will be computed.
1182
1183 Returns
1184 -------
1185 w : (..., M,) ndarray
1186 The eigenvalues, each repeated according to its multiplicity.
1187 They are not necessarily ordered, nor are they necessarily
1188 real for real matrices.
1189
1190 Raises
1191 ------
1192 LinAlgError
1193 If the eigenvalue computation does not converge.
1194
1195 See Also
1196 --------
1197 eig : eigenvalues and right eigenvectors of general arrays
1198 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1199 (conjugate symmetric) arrays.
1200 eigh : eigenvalues and eigenvectors of real symmetric or complex