linalg.py 87 KB

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