test_linalg.py 73 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089
  1. """ Test functions for linalg module
  2. """
  3. import os
  4. import sys
  5. import itertools
  6. import traceback
  7. import textwrap
  8. import subprocess
  9. import pytest
  10. import numpy as np
  11. from numpy import array, single, double, csingle, cdouble, dot, identity, matmul
  12. from numpy import multiply, atleast_2d, inf, asarray
  13. from numpy import linalg
  14. from numpy.linalg import matrix_power, norm, matrix_rank, multi_dot, LinAlgError
  15. from numpy.linalg.linalg import _multi_dot_matrix_chain_order
  16. from numpy.testing import (
  17. assert_, assert_equal, assert_raises, assert_array_equal,
  18. assert_almost_equal, assert_allclose, suppress_warnings,
  19. assert_raises_regex, HAS_LAPACK64,
  20. )
  21. from numpy.testing._private.utils import requires_memory
  22. def consistent_subclass(out, in_):
  23. # For ndarray subclass input, our output should have the same subclass
  24. # (non-ndarray input gets converted to ndarray).
  25. return type(out) is (type(in_) if isinstance(in_, np.ndarray)
  26. else np.ndarray)
  27. old_assert_almost_equal = assert_almost_equal
  28. def assert_almost_equal(a, b, single_decimal=6, double_decimal=12, **kw):
  29. if asarray(a).dtype.type in (single, csingle):
  30. decimal = single_decimal
  31. else:
  32. decimal = double_decimal
  33. old_assert_almost_equal(a, b, decimal=decimal, **kw)
  34. def get_real_dtype(dtype):
  35. return {single: single, double: double,
  36. csingle: single, cdouble: double}[dtype]
  37. def get_complex_dtype(dtype):
  38. return {single: csingle, double: cdouble,
  39. csingle: csingle, cdouble: cdouble}[dtype]
  40. def get_rtol(dtype):
  41. # Choose a safe rtol
  42. if dtype in (single, csingle):
  43. return 1e-5
  44. else:
  45. return 1e-11
  46. # used to categorize tests
  47. all_tags = {
  48. 'square', 'nonsquare', 'hermitian', # mutually exclusive
  49. 'generalized', 'size-0', 'strided' # optional additions
  50. }
  51. class LinalgCase:
  52. def __init__(self, name, a, b, tags=set()):
  53. """
  54. A bundle of arguments to be passed to a test case, with an identifying
  55. name, the operands a and b, and a set of tags to filter the tests
  56. """
  57. assert_(isinstance(name, str))
  58. self.name = name
  59. self.a = a
  60. self.b = b
  61. self.tags = frozenset(tags) # prevent shared tags
  62. def check(self, do):
  63. """
  64. Run the function `do` on this test case, expanding arguments
  65. """
  66. do(self.a, self.b, tags=self.tags)
  67. def __repr__(self):
  68. return f'<LinalgCase: {self.name}>'
  69. def apply_tag(tag, cases):
  70. """
  71. Add the given tag (a string) to each of the cases (a list of LinalgCase
  72. objects)
  73. """
  74. assert tag in all_tags, "Invalid tag"
  75. for case in cases:
  76. case.tags = case.tags | {tag}
  77. return cases
  78. #
  79. # Base test cases
  80. #
  81. np.random.seed(1234)
  82. CASES = []
  83. # square test cases
  84. CASES += apply_tag('square', [
  85. LinalgCase("single",
  86. array([[1., 2.], [3., 4.]], dtype=single),
  87. array([2., 1.], dtype=single)),
  88. LinalgCase("double",
  89. array([[1., 2.], [3., 4.]], dtype=double),
  90. array([2., 1.], dtype=double)),
  91. LinalgCase("double_2",
  92. array([[1., 2.], [3., 4.]], dtype=double),
  93. array([[2., 1., 4.], [3., 4., 6.]], dtype=double)),
  94. LinalgCase("csingle",
  95. array([[1. + 2j, 2 + 3j], [3 + 4j, 4 + 5j]], dtype=csingle),
  96. array([2. + 1j, 1. + 2j], dtype=csingle)),
  97. LinalgCase("cdouble",
  98. array([[1. + 2j, 2 + 3j], [3 + 4j, 4 + 5j]], dtype=cdouble),
  99. array([2. + 1j, 1. + 2j], dtype=cdouble)),
  100. LinalgCase("cdouble_2",
  101. array([[1. + 2j, 2 + 3j], [3 + 4j, 4 + 5j]], dtype=cdouble),
  102. array([[2. + 1j, 1. + 2j, 1 + 3j], [1 - 2j, 1 - 3j, 1 - 6j]], dtype=cdouble)),
  103. LinalgCase("0x0",
  104. np.empty((0, 0), dtype=double),
  105. np.empty((0,), dtype=double),
  106. tags={'size-0'}),
  107. LinalgCase("8x8",
  108. np.random.rand(8, 8),
  109. np.random.rand(8)),
  110. LinalgCase("1x1",
  111. np.random.rand(1, 1),
  112. np.random.rand(1)),
  113. LinalgCase("nonarray",
  114. [[1, 2], [3, 4]],
  115. [2, 1]),
  116. ])
  117. # non-square test-cases
  118. CASES += apply_tag('nonsquare', [
  119. LinalgCase("single_nsq_1",
  120. array([[1., 2., 3.], [3., 4., 6.]], dtype=single),
  121. array([2., 1.], dtype=single)),
  122. LinalgCase("single_nsq_2",
  123. array([[1., 2.], [3., 4.], [5., 6.]], dtype=single),
  124. array([2., 1., 3.], dtype=single)),
  125. LinalgCase("double_nsq_1",
  126. array([[1., 2., 3.], [3., 4., 6.]], dtype=double),
  127. array([2., 1.], dtype=double)),
  128. LinalgCase("double_nsq_2",
  129. array([[1., 2.], [3., 4.], [5., 6.]], dtype=double),
  130. array([2., 1., 3.], dtype=double)),
  131. LinalgCase("csingle_nsq_1",
  132. array(
  133. [[1. + 1j, 2. + 2j, 3. - 3j], [3. - 5j, 4. + 9j, 6. + 2j]], dtype=csingle),
  134. array([2. + 1j, 1. + 2j], dtype=csingle)),
  135. LinalgCase("csingle_nsq_2",
  136. array(
  137. [[1. + 1j, 2. + 2j], [3. - 3j, 4. - 9j], [5. - 4j, 6. + 8j]], dtype=csingle),
  138. array([2. + 1j, 1. + 2j, 3. - 3j], dtype=csingle)),
  139. LinalgCase("cdouble_nsq_1",
  140. array(
  141. [[1. + 1j, 2. + 2j, 3. - 3j], [3. - 5j, 4. + 9j, 6. + 2j]], dtype=cdouble),
  142. array([2. + 1j, 1. + 2j], dtype=cdouble)),
  143. LinalgCase("cdouble_nsq_2",
  144. array(
  145. [[1. + 1j, 2. + 2j], [3. - 3j, 4. - 9j], [5. - 4j, 6. + 8j]], dtype=cdouble),
  146. array([2. + 1j, 1. + 2j, 3. - 3j], dtype=cdouble)),
  147. LinalgCase("cdouble_nsq_1_2",
  148. array(
  149. [[1. + 1j, 2. + 2j, 3. - 3j], [3. - 5j, 4. + 9j, 6. + 2j]], dtype=cdouble),
  150. array([[2. + 1j, 1. + 2j], [1 - 1j, 2 - 2j]], dtype=cdouble)),
  151. LinalgCase("cdouble_nsq_2_2",
  152. array(
  153. [[1. + 1j, 2. + 2j], [3. - 3j, 4. - 9j], [5. - 4j, 6. + 8j]], dtype=cdouble),
  154. array([[2. + 1j, 1. + 2j], [1 - 1j, 2 - 2j], [1 - 1j, 2 - 2j]], dtype=cdouble)),
  155. LinalgCase("8x11",
  156. np.random.rand(8, 11),
  157. np.random.rand(8)),
  158. LinalgCase("1x5",
  159. np.random.rand(1, 5),
  160. np.random.rand(1)),
  161. LinalgCase("5x1",
  162. np.random.rand(5, 1),
  163. np.random.rand(5)),
  164. LinalgCase("0x4",
  165. np.random.rand(0, 4),
  166. np.random.rand(0),
  167. tags={'size-0'}),
  168. LinalgCase("4x0",
  169. np.random.rand(4, 0),
  170. np.random.rand(4),
  171. tags={'size-0'}),
  172. ])
  173. # hermitian test-cases
  174. CASES += apply_tag('hermitian', [
  175. LinalgCase("hsingle",
  176. array([[1., 2.], [2., 1.]], dtype=single),
  177. None),
  178. LinalgCase("hdouble",
  179. array([[1., 2.], [2., 1.]], dtype=double),
  180. None),
  181. LinalgCase("hcsingle",
  182. array([[1., 2 + 3j], [2 - 3j, 1]], dtype=csingle),
  183. None),
  184. LinalgCase("hcdouble",
  185. array([[1., 2 + 3j], [2 - 3j, 1]], dtype=cdouble),
  186. None),
  187. LinalgCase("hempty",
  188. np.empty((0, 0), dtype=double),
  189. None,
  190. tags={'size-0'}),
  191. LinalgCase("hnonarray",
  192. [[1, 2], [2, 1]],
  193. None),
  194. LinalgCase("matrix_b_only",
  195. array([[1., 2.], [2., 1.]]),
  196. None),
  197. LinalgCase("hmatrix_1x1",
  198. np.random.rand(1, 1),
  199. None),
  200. ])
  201. #
  202. # Gufunc test cases
  203. #
  204. def _make_generalized_cases():
  205. new_cases = []
  206. for case in CASES:
  207. if not isinstance(case.a, np.ndarray):
  208. continue
  209. a = np.array([case.a, 2 * case.a, 3 * case.a])
  210. if case.b is None:
  211. b = None
  212. else:
  213. b = np.array([case.b, 7 * case.b, 6 * case.b])
  214. new_case = LinalgCase(case.name + "_tile3", a, b,
  215. tags=case.tags | {'generalized'})
  216. new_cases.append(new_case)
  217. a = np.array([case.a] * 2 * 3).reshape((3, 2) + case.a.shape)
  218. if case.b is None:
  219. b = None
  220. else:
  221. b = np.array([case.b] * 2 * 3).reshape((3, 2) + case.b.shape)
  222. new_case = LinalgCase(case.name + "_tile213", a, b,
  223. tags=case.tags | {'generalized'})
  224. new_cases.append(new_case)
  225. return new_cases
  226. CASES += _make_generalized_cases()
  227. #
  228. # Generate stride combination variations of the above
  229. #
  230. def _stride_comb_iter(x):
  231. """
  232. Generate cartesian product of strides for all axes
  233. """
  234. if not isinstance(x, np.ndarray):
  235. yield x, "nop"
  236. return
  237. stride_set = [(1,)] * x.ndim
  238. stride_set[-1] = (1, 3, -4)
  239. if x.ndim > 1:
  240. stride_set[-2] = (1, 3, -4)
  241. if x.ndim > 2:
  242. stride_set[-3] = (1, -4)
  243. for repeats in itertools.product(*tuple(stride_set)):
  244. new_shape = [abs(a * b) for a, b in zip(x.shape, repeats)]
  245. slices = tuple([slice(None, None, repeat) for repeat in repeats])
  246. # new array with different strides, but same data
  247. xi = np.empty(new_shape, dtype=x.dtype)
  248. xi.view(np.uint32).fill(0xdeadbeef)
  249. xi = xi[slices]
  250. xi[...] = x
  251. xi = xi.view(x.__class__)
  252. assert_(np.all(xi == x))
  253. yield xi, "stride_" + "_".join(["%+d" % j for j in repeats])
  254. # generate also zero strides if possible
  255. if x.ndim >= 1 and x.shape[-1] == 1:
  256. s = list(x.strides)
  257. s[-1] = 0
  258. xi = np.lib.stride_tricks.as_strided(x, strides=s)
  259. yield xi, "stride_xxx_0"
  260. if x.ndim >= 2 and x.shape[-2] == 1:
  261. s = list(x.strides)
  262. s[-2] = 0
  263. xi = np.lib.stride_tricks.as_strided(x, strides=s)
  264. yield xi, "stride_xxx_0_x"
  265. if x.ndim >= 2 and x.shape[:-2] == (1, 1):
  266. s = list(x.strides)
  267. s[-1] = 0
  268. s[-2] = 0
  269. xi = np.lib.stride_tricks.as_strided(x, strides=s)
  270. yield xi, "stride_xxx_0_0"
  271. def _make_strided_cases():
  272. new_cases = []
  273. for case in CASES:
  274. for a, a_label in _stride_comb_iter(case.a):
  275. for b, b_label in _stride_comb_iter(case.b):
  276. new_case = LinalgCase(case.name + "_" + a_label + "_" + b_label, a, b,
  277. tags=case.tags | {'strided'})
  278. new_cases.append(new_case)
  279. return new_cases
  280. CASES += _make_strided_cases()
  281. #
  282. # Test different routines against the above cases
  283. #
  284. class LinalgTestCase:
  285. TEST_CASES = CASES
  286. def check_cases(self, require=set(), exclude=set()):
  287. """
  288. Run func on each of the cases with all of the tags in require, and none
  289. of the tags in exclude
  290. """
  291. for case in self.TEST_CASES:
  292. # filter by require and exclude
  293. if case.tags & require != require:
  294. continue
  295. if case.tags & exclude:
  296. continue
  297. try:
  298. case.check(self.do)
  299. except Exception:
  300. msg = f'In test case: {case!r}\n\n'
  301. msg += traceback.format_exc()
  302. raise AssertionError(msg)
  303. class LinalgSquareTestCase(LinalgTestCase):
  304. def test_sq_cases(self):
  305. self.check_cases(require={'square'},
  306. exclude={'generalized', 'size-0'})
  307. def test_empty_sq_cases(self):
  308. self.check_cases(require={'square', 'size-0'},
  309. exclude={'generalized'})
  310. class LinalgNonsquareTestCase(LinalgTestCase):
  311. def test_nonsq_cases(self):
  312. self.check_cases(require={'nonsquare'},
  313. exclude={'generalized', 'size-0'})
  314. def test_empty_nonsq_cases(self):
  315. self.check_cases(require={'nonsquare', 'size-0'},
  316. exclude={'generalized'})
  317. class HermitianTestCase(LinalgTestCase):
  318. def test_herm_cases(self):
  319. self.check_cases(require={'hermitian'},
  320. exclude={'generalized', 'size-0'})
  321. def test_empty_herm_cases(self):
  322. self.check_cases(require={'hermitian', 'size-0'},
  323. exclude={'generalized'})
  324. class LinalgGeneralizedSquareTestCase(LinalgTestCase):
  325. @pytest.mark.slow
  326. def test_generalized_sq_cases(self):
  327. self.check_cases(require={'generalized', 'square'},
  328. exclude={'size-0'})
  329. @pytest.mark.slow
  330. def test_generalized_empty_sq_cases(self):
  331. self.check_cases(require={'generalized', 'square', 'size-0'})
  332. class LinalgGeneralizedNonsquareTestCase(LinalgTestCase):
  333. @pytest.mark.slow
  334. def test_generalized_nonsq_cases(self):
  335. self.check_cases(require={'generalized', 'nonsquare'},
  336. exclude={'size-0'})
  337. @pytest.mark.slow
  338. def test_generalized_empty_nonsq_cases(self):
  339. self.check_cases(require={'generalized', 'nonsquare', 'size-0'})
  340. class HermitianGeneralizedTestCase(LinalgTestCase):
  341. @pytest.mark.slow
  342. def test_generalized_herm_cases(self):
  343. self.check_cases(require={'generalized', 'hermitian'},
  344. exclude={'size-0'})
  345. @pytest.mark.slow
  346. def test_generalized_empty_herm_cases(self):
  347. self.check_cases(require={'generalized', 'hermitian', 'size-0'},
  348. exclude={'none'})
  349. def dot_generalized(a, b):
  350. a = asarray(a)
  351. if a.ndim >= 3:
  352. if a.ndim == b.ndim:
  353. # matrix x matrix
  354. new_shape = a.shape[:-1] + b.shape[-1:]
  355. elif a.ndim == b.ndim + 1:
  356. # matrix x vector
  357. new_shape = a.shape[:-1]
  358. else:
  359. raise ValueError("Not implemented...")
  360. r = np.empty(new_shape, dtype=np.common_type(a, b))
  361. for c in itertools.product(*map(range, a.shape[:-2])):
  362. r[c] = dot(a[c], b[c])
  363. return r
  364. else:
  365. return dot(a, b)
  366. def identity_like_generalized(a):
  367. a = asarray(a)
  368. if a.ndim >= 3:
  369. r = np.empty(a.shape, dtype=a.dtype)
  370. r[...] = identity(a.shape[-2])
  371. return r
  372. else:
  373. return identity(a.shape[0])
  374. class SolveCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
  375. # kept apart from TestSolve for use for testing with matrices.
  376. def do(self, a, b, tags):
  377. x = linalg.solve(a, b)
  378. assert_almost_equal(b, dot_generalized(a, x))
  379. assert_(consistent_subclass(x, b))
  380. class TestSolve(SolveCases):
  381. @pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
  382. def test_types(self, dtype):
  383. x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
  384. assert_equal(linalg.solve(x, x).dtype, dtype)
  385. def test_0_size(self):
  386. class ArraySubclass(np.ndarray):
  387. pass
  388. # Test system of 0x0 matrices
  389. a = np.arange(8).reshape(2, 2, 2)
  390. b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)
  391. expected = linalg.solve(a, b)[:, 0:0, :]
  392. result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, :])
  393. assert_array_equal(result, expected)
  394. assert_(isinstance(result, ArraySubclass))
  395. # Test errors for non-square and only b's dimension being 0
  396. assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
  397. assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])
  398. # Test broadcasting error
  399. b = np.arange(6).reshape(1, 3, 2) # broadcasting error
  400. assert_raises(ValueError, linalg.solve, a, b)
  401. assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
  402. # Test zero "single equations" with 0x0 matrices.
  403. b = np.arange(2).reshape(1, 2).view(ArraySubclass)
  404. expected = linalg.solve(a, b)[:, 0:0]
  405. result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
  406. assert_array_equal(result, expected)
  407. assert_(isinstance(result, ArraySubclass))
  408. b = np.arange(3).reshape(1, 3)
  409. assert_raises(ValueError, linalg.solve, a, b)
  410. assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
  411. assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b)
  412. def test_0_size_k(self):
  413. # test zero multiple equation (K=0) case.
  414. class ArraySubclass(np.ndarray):
  415. pass
  416. a = np.arange(4).reshape(1, 2, 2)
  417. b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass)
  418. expected = linalg.solve(a, b)[:, :, 0:0]
  419. result = linalg.solve(a, b[:, :, 0:0])
  420. assert_array_equal(result, expected)
  421. assert_(isinstance(result, ArraySubclass))
  422. # test both zero.
  423. expected = linalg.solve(a, b)[:, 0:0, 0:0]
  424. result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, 0:0])
  425. assert_array_equal(result, expected)
  426. assert_(isinstance(result, ArraySubclass))
  427. class InvCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
  428. def do(self, a, b, tags):
  429. a_inv = linalg.inv(a)
  430. assert_almost_equal(dot_generalized(a, a_inv),
  431. identity_like_generalized(a))
  432. assert_(consistent_subclass(a_inv, a))
  433. class TestInv(InvCases):
  434. @pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
  435. def test_types(self, dtype):
  436. x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
  437. assert_equal(linalg.inv(x).dtype, dtype)
  438. def test_0_size(self):
  439. # Check that all kinds of 0-sized arrays work
  440. class ArraySubclass(np.ndarray):
  441. pass
  442. a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
  443. res = linalg.inv(a)
  444. assert_(res.dtype.type is np.float64)
  445. assert_equal(a.shape, res.shape)
  446. assert_(isinstance(res, ArraySubclass))
  447. a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
  448. res = linalg.inv(a)
  449. assert_(res.dtype.type is np.complex64)
  450. assert_equal(a.shape, res.shape)
  451. assert_(isinstance(res, ArraySubclass))
  452. class EigvalsCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
  453. def do(self, a, b, tags):
  454. ev = linalg.eigvals(a)
  455. evalues, evectors = linalg.eig(a)
  456. assert_almost_equal(ev, evalues)
  457. class TestEigvals(EigvalsCases):
  458. @pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
  459. def test_types(self, dtype):
  460. x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
  461. assert_equal(linalg.eigvals(x).dtype, dtype)
  462. x = np.array([[1, 0.5], [-1, 1]], dtype=dtype)
  463. assert_equal(linalg.eigvals(x).dtype, get_complex_dtype(dtype))
  464. def test_0_size(self):
  465. # Check that all kinds of 0-sized arrays work
  466. class ArraySubclass(np.ndarray):
  467. pass
  468. a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
  469. res = linalg.eigvals(a)
  470. assert_(res.dtype.type is np.float64)
  471. assert_equal((0, 1), res.shape)
  472. # This is just for documentation, it might make sense to change:
  473. assert_(isinstance(res, np.ndarray))
  474. a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
  475. res = linalg.eigvals(a)
  476. assert_(res.dtype.type is np.complex64)
  477. assert_equal((0,), res.shape)
  478. # This is just for documentation, it might make sense to change:
  479. assert_(isinstance(res, np.ndarray))
  480. class EigCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
  481. def do(self, a, b, tags):
  482. evalues, evectors = linalg.eig(a)
  483. assert_allclose(dot_generalized(a, evectors),
  484. np.asarray(evectors) * np.asarray(evalues)[..., None, :],
  485. rtol=get_rtol(evalues.dtype))
  486. assert_(consistent_subclass(evectors, a))
  487. class TestEig(EigCases):
  488. @pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
  489. def test_types(self, dtype):
  490. x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
  491. w, v = np.linalg.eig(x)
  492. assert_equal(w.dtype, dtype)
  493. assert_equal(v.dtype, dtype)
  494. x = np.array([[1, 0.5], [-1, 1]], dtype=dtype)
  495. w, v = np.linalg.eig(x)
  496. assert_equal(w.dtype, get_complex_dtype(dtype))
  497. assert_equal(v.dtype, get_complex_dtype(dtype))
  498. def test_0_size(self):
  499. # Check that all kinds of 0-sized arrays work
  500. class ArraySubclass(np.ndarray):
  501. pass
  502. a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
  503. res, res_v = linalg.eig(a)
  504. assert_(res_v.dtype.type is np.float64)
  505. assert_(res.dtype.type is np.float64)
  506. assert_equal(a.shape, res_v.shape)
  507. assert_equal((0, 1), res.shape)
  508. # This is just for documentation, it might make sense to change:
  509. assert_(isinstance(a, np.ndarray))
  510. a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
  511. res, res_v = linalg.eig(a)
  512. assert_(res_v.dtype.type is np.complex64)
  513. assert_(res.dtype.type is np.complex64)
  514. assert_equal(a.shape, res_v.shape)
  515. assert_equal((0,), res.shape)
  516. # This is just for documentation, it might make sense to change:
  517. assert_(isinstance(a, np.ndarray))
  518. class SVDBaseTests:
  519. hermitian = False
  520. @pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
  521. def test_types(self, dtype):
  522. x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
  523. u, s, vh = linalg.svd(x)
  524. assert_equal(u.dtype, dtype)
  525. assert_equal(s.dtype, get_real_dtype(dtype))
  526. assert_equal(vh.dtype, dtype)
  527. s = linalg.svd(x, compute_uv=False, hermitian=self.hermitian)
  528. assert_equal(s.dtype, get_real_dtype(dtype))
  529. class SVDCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
  530. def do(self, a, b, tags):
  531. u, s, vt = linalg.svd(a, False)
  532. assert_allclose(a, dot_generalized(np.asarray(u) * np.asarray(s)[..., None, :],
  533. np.asarray(vt)),
  534. rtol=get_rtol(u.dtype))
  535. assert_(consistent_subclass(u, a))
  536. assert_(consistent_subclass(vt, a))
  537. class TestSVD(SVDCases, SVDBaseTests):
  538. def test_empty_identity(self):
  539. """ Empty input should put an identity matrix in u or vh """
  540. x = np.empty((4, 0))
  541. u, s, vh = linalg.svd(x, compute_uv=True, hermitian=self.hermitian)
  542. assert_equal(u.shape, (4, 4))
  543. assert_equal(vh.shape, (0, 0))
  544. assert_equal(u, np.eye(4))
  545. x = np.empty((0, 4))
  546. u, s, vh = linalg.svd(x, compute_uv=True, hermitian=self.hermitian)
  547. assert_equal(u.shape, (0, 0))
  548. assert_equal(vh.shape, (4, 4))
  549. assert_equal(vh, np.eye(4))
  550. class SVDHermitianCases(HermitianTestCase, HermitianGeneralizedTestCase):
  551. def do(self, a, b, tags):
  552. u, s, vt = linalg.svd(a, False, hermitian=True)
  553. assert_allclose(a, dot_generalized(np.asarray(u) * np.asarray(s)[..., None, :],
  554. np.asarray(vt)),
  555. rtol=get_rtol(u.dtype))
  556. def hermitian(mat):
  557. axes = list(range(mat.ndim))
  558. axes[-1], axes[-2] = axes[-2], axes[-1]
  559. return np.conj(np.transpose(mat, axes=axes))
  560. assert_almost_equal(np.matmul(u, hermitian(u)), np.broadcast_to(np.eye(u.shape[-1]), u.shape))
  561. assert_almost_equal(np.matmul(vt, hermitian(vt)), np.broadcast_to(np.eye(vt.shape[-1]), vt.shape))
  562. assert_equal(np.sort(s)[..., ::-1], s)
  563. assert_(consistent_subclass(u, a))
  564. assert_(consistent_subclass(vt, a))
  565. class TestSVDHermitian(SVDHermitianCases, SVDBaseTests):
  566. hermitian = True
  567. class CondCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
  568. # cond(x, p) for p in (None, 2, -2)
  569. def do(self, a, b, tags):
  570. c = asarray(a) # a might be a matrix
  571. if 'size-0' in tags:
  572. assert_raises(LinAlgError, linalg.cond, c)
  573. return
  574. # +-2 norms
  575. s = linalg.svd(c, compute_uv=False)
  576. assert_almost_equal(
  577. linalg.cond(a), s[..., 0] / s[..., -1],
  578. single_decimal=5, double_decimal=11)
  579. assert_almost_equal(
  580. linalg.cond(a, 2), s[..., 0] / s[..., -1],
  581. single_decimal=5, double_decimal=11)
  582. assert_almost_equal(
  583. linalg.cond(a, -2), s[..., -1] / s[..., 0],
  584. single_decimal=5, double_decimal=11)
  585. # Other norms
  586. cinv = np.linalg.inv(c)
  587. assert_almost_equal(
  588. linalg.cond(a, 1),
  589. abs(c).sum(-2).max(-1) * abs(cinv).sum(-2).max(-1),
  590. single_decimal=5, double_decimal=11)
  591. assert_almost_equal(
  592. linalg.cond(a, -1),
  593. abs(c).sum(-2).min(-1) * abs(cinv).sum(-2).min(-1),
  594. single_decimal=5, double_decimal=11)
  595. assert_almost_equal(
  596. linalg.cond(a, np.inf),
  597. abs(c).sum(-1).max(-1) * abs(cinv).sum(-1).max(-1),
  598. single_decimal=5, double_decimal=11)
  599. assert_almost_equal(
  600. linalg.cond(a, -np.inf),
  601. abs(c).sum(-1).min(-1) * abs(cinv).sum(-1).min(-1),
  602. single_decimal=5, double_decimal=11)
  603. assert_almost_equal(
  604. linalg.cond(a, 'fro'),
  605. np.sqrt((abs(c)**2).sum(-1).sum(-1)
  606. * (abs(cinv)**2).sum(-1).sum(-1)),
  607. single_decimal=5, double_decimal=11)
  608. class TestCond(CondCases):
  609. def test_basic_nonsvd(self):
  610. # Smoketest the non-svd norms
  611. A = array([[1., 0, 1], [0, -2., 0], [0, 0, 3.]])
  612. assert_almost_equal(linalg.cond(A, inf), 4)
  613. assert_almost_equal(linalg.cond(A, -inf), 2/3)
  614. assert_almost_equal(linalg.cond(A, 1), 4)
  615. assert_almost_equal(linalg.cond(A, -1), 0.5)
  616. assert_almost_equal(linalg.cond(A, 'fro'), np.sqrt(265 / 12))
  617. def test_singular(self):
  618. # Singular matrices have infinite condition number for
  619. # positive norms, and negative norms shouldn't raise
  620. # exceptions
  621. As = [np.zeros((2, 2)), np.ones((2, 2))]
  622. p_pos = [None, 1, 2, 'fro']
  623. p_neg = [-1, -2]
  624. for A, p in itertools.product(As, p_pos):
  625. # Inversion may not hit exact infinity, so just check the
  626. # number is large
  627. assert_(linalg.cond(A, p) > 1e15)
  628. for A, p in itertools.product(As, p_neg):
  629. linalg.cond(A, p)
  630. def test_nan(self):
  631. # nans should be passed through, not converted to infs
  632. ps = [None, 1, -1, 2, -2, 'fro']
  633. p_pos = [None, 1, 2, 'fro']
  634. A = np.ones((2, 2))
  635. A[0,1] = np.nan
  636. for p in ps:
  637. c = linalg.cond(A, p)
  638. assert_(isinstance(c, np.float_))
  639. assert_(np.isnan(c))
  640. A = np.ones((3, 2, 2))
  641. A[1,0,1] = np.nan
  642. for p in ps:
  643. c = linalg.cond(A, p)
  644. assert_(np.isnan(c[1]))
  645. if p in p_pos:
  646. assert_(c[0] > 1e15)
  647. assert_(c[2] > 1e15)
  648. else:
  649. assert_(not np.isnan(c[0]))
  650. assert_(not np.isnan(c[2]))
  651. def test_stacked_singular(self):
  652. # Check behavior when only some of the stacked matrices are
  653. # singular
  654. np.random.seed(1234)
  655. A = np.random.rand(2, 2, 2, 2)
  656. A[0,0] = 0
  657. A[1,1] = 0
  658. for p in (None, 1, 2, 'fro', -1, -2):
  659. c = linalg.cond(A, p)
  660. assert_equal(c[0,0], np.inf)
  661. assert_equal(c[1,1], np.inf)
  662. assert_(np.isfinite(c[0,1]))
  663. assert_(np.isfinite(c[1,0]))
  664. class PinvCases(LinalgSquareTestCase,
  665. LinalgNonsquareTestCase,
  666. LinalgGeneralizedSquareTestCase,
  667. LinalgGeneralizedNonsquareTestCase):
  668. def do(self, a, b, tags):
  669. a_ginv = linalg.pinv(a)
  670. # `a @ a_ginv == I` does not hold if a is singular
  671. dot = dot_generalized
  672. assert_almost_equal(dot(dot(a, a_ginv), a), a, single_decimal=5, double_decimal=11)
  673. assert_(consistent_subclass(a_ginv, a))
  674. class TestPinv(PinvCases):
  675. pass
  676. class PinvHermitianCases(HermitianTestCase, HermitianGeneralizedTestCase):
  677. def do(self, a, b, tags):
  678. a_ginv = linalg.pinv(a, hermitian=True)
  679. # `a @ a_ginv == I` does not hold if a is singular
  680. dot = dot_generalized
  681. assert_almost_equal(dot(dot(a, a_ginv), a), a, single_decimal=5, double_decimal=11)
  682. assert_(consistent_subclass(a_ginv, a))
  683. class TestPinvHermitian(PinvHermitianCases):
  684. pass
  685. class DetCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
  686. def do(self, a, b, tags):
  687. d = linalg.det(a)
  688. (s, ld) = linalg.slogdet(a)
  689. if asarray(a).dtype.type in (single, double):
  690. ad = asarray(a).astype(double)
  691. else:
  692. ad = asarray(a).astype(cdouble)
  693. ev = linalg.eigvals(ad)
  694. assert_almost_equal(d, multiply.reduce(ev, axis=-1))
  695. assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1))
  696. s = np.atleast_1d(s)
  697. ld = np.atleast_1d(ld)
  698. m = (s != 0)
  699. assert_almost_equal(np.abs(s[m]), 1)
  700. assert_equal(ld[~m], -inf)
  701. class TestDet(DetCases):
  702. def test_zero(self):
  703. assert_equal(linalg.det([[0.0]]), 0.0)
  704. assert_equal(type(linalg.det([[0.0]])), double)
  705. assert_equal(linalg.det([[0.0j]]), 0.0)
  706. assert_equal(type(linalg.det([[0.0j]])), cdouble)
  707. assert_equal(linalg.slogdet([[0.0]]), (0.0, -inf))
  708. assert_equal(type(linalg.slogdet([[0.0]])[0]), double)
  709. assert_equal(type(linalg.slogdet([[0.0]])[1]), double)
  710. assert_equal(linalg.slogdet([[0.0j]]), (0.0j, -inf))
  711. assert_equal(type(linalg.slogdet([[0.0j]])[0]), cdouble)
  712. assert_equal(type(linalg.slogdet([[0.0j]])[1]), double)
  713. @pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
  714. def test_types(self, dtype):
  715. x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
  716. assert_equal(np.linalg.det(x).dtype, dtype)
  717. ph, s = np.linalg.slogdet(x)
  718. assert_equal(s.dtype, get_real_dtype(dtype))
  719. assert_equal(ph.dtype, dtype)
  720. def test_0_size(self):
  721. a = np.zeros((0, 0), dtype=np.complex64)
  722. res = linalg.det(a)
  723. assert_equal(res, 1.)
  724. assert_(res.dtype.type is np.complex64)
  725. res = linalg.slogdet(a)
  726. assert_equal(res, (1, 0))
  727. assert_(res[0].dtype.type is np.complex64)
  728. assert_(res[1].dtype.type is np.float32)
  729. a = np.zeros((0, 0), dtype=np.float64)
  730. res = linalg.det(a)
  731. assert_equal(res, 1.)
  732. assert_(res.dtype.type is np.float64)
  733. res = linalg.slogdet(a)
  734. assert_equal(res, (1, 0))
  735. assert_(res[0].dtype.type is np.float64)
  736. assert_(res[1].dtype.type is np.float64)
  737. class LstsqCases(LinalgSquareTestCase, LinalgNonsquareTestCase):
  738. def do(self, a, b, tags):
  739. arr = np.asarray(a)
  740. m, n = arr.shape
  741. u, s, vt = linalg.svd(a, False)
  742. x, residuals, rank, sv = linalg.lstsq(a, b, rcond=-1)
  743. if m == 0:
  744. assert_((x == 0).all())
  745. if m <= n:
  746. assert_almost_equal(b, dot(a, x))
  747. assert_equal(rank, m)
  748. else:
  749. assert_equal(rank, n)
  750. assert_almost_equal(sv, sv.__array_wrap__(s))
  751. if rank == n and m > n:
  752. expect_resids = (
  753. np.asarray(abs(np.dot(a, x) - b)) ** 2).sum(axis=0)
  754. expect_resids = np.asarray(expect_resids)
  755. if np.asarray(b).ndim == 1:
  756. expect_resids.shape = (1,)
  757. assert_equal(residuals.shape, expect_resids.shape)
  758. else:
  759. expect_resids = np.array([]).view(type(x))
  760. assert_almost_equal(residuals, expect_resids)
  761. assert_(np.issubdtype(residuals.dtype, np.floating))
  762. assert_(consistent_subclass(x, b))
  763. assert_(consistent_subclass(residuals, b))
  764. class TestLstsq(LstsqCases):
  765. def test_future_rcond(self):
  766. a = np.array([[0., 1., 0., 1., 2., 0.],
  767. [0., 2., 0., 0., 1., 0.],
  768. [1., 0., 1., 0., 0., 4.],
  769. [0., 0., 0., 2., 3., 0.]]).T
  770. b = np.array([1, 0, 0, 0, 0, 0])
  771. with suppress_warnings() as sup:
  772. w = sup.record(FutureWarning, "`rcond` parameter will change")
  773. x, residuals, rank, s = linalg.lstsq(a, b)
  774. assert_(rank == 4)
  775. x, residuals, rank, s = linalg.lstsq(a, b, rcond=-1)
  776. assert_(rank == 4)
  777. x, residuals, rank, s = linalg.lstsq(a, b, rcond=None)
  778. assert_(rank == 3)
  779. # Warning should be raised exactly once (first command)
  780. assert_(len(w) == 1)
  781. @pytest.mark.parametrize(["m", "n", "n_rhs"], [
  782. (4, 2, 2),
  783. (0, 4, 1),
  784. (0, 4, 2),
  785. (4, 0, 1),
  786. (4, 0, 2),
  787. (4, 2, 0),
  788. (0, 0, 0)
  789. ])
  790. def test_empty_a_b(self, m, n, n_rhs):
  791. a = np.arange(m * n).reshape(m, n)
  792. b = np.ones((m, n_rhs))
  793. x, residuals, rank, s = linalg.lstsq(a, b, rcond=None)
  794. if m == 0:
  795. assert_((x == 0).all())
  796. assert_equal(x.shape, (n, n_rhs))
  797. assert_equal(residuals.shape, ((n_rhs,) if m > n else (0,)))
  798. if m > n and n_rhs > 0:
  799. # residuals are exactly the squared norms of b's columns
  800. r = b - np.dot(a, x)
  801. assert_almost_equal(residuals, (r * r).sum(axis=-2))
  802. assert_equal(rank, min(m, n))
  803. assert_equal(s.shape, (min(m, n),))
  804. def test_incompatible_dims(self):
  805. # use modified version of docstring example
  806. x = np.array([0, 1, 2, 3])
  807. y = np.array([-1, 0.2, 0.9, 2.1, 3.3])
  808. A = np.vstack([x, np.ones(len(x))]).T
  809. with assert_raises_regex(LinAlgError, "Incompatible dimensions"):
  810. linalg.lstsq(A, y, rcond=None)
  811. @pytest.mark.parametrize('dt', [np.dtype(c) for c in '?bBhHiIqQefdgFDGO'])
  812. class TestMatrixPower:
  813. rshft_0 = np.eye(4)
  814. rshft_1 = rshft_0[[3, 0, 1, 2]]
  815. rshft_2 = rshft_0[[2, 3, 0, 1]]
  816. rshft_3 = rshft_0[[1, 2, 3, 0]]
  817. rshft_all = [rshft_0, rshft_1, rshft_2, rshft_3]
  818. noninv = array([[1, 0], [0, 0]])
  819. stacked = np.block([[[rshft_0]]]*2)
  820. #FIXME the 'e' dtype might work in future
  821. dtnoinv = [object, np.dtype('e'), np.dtype('g'), np.dtype('G')]
  822. def test_large_power(self, dt):
  823. rshft = self.rshft_1.astype(dt)
  824. assert_equal(
  825. matrix_power(rshft, 2**100 + 2**10 + 2**5 + 0), self.rshft_0)
  826. assert_equal(
  827. matrix_power(rshft, 2**100 + 2**10 + 2**5 + 1), self.rshft_1)
  828. assert_equal(
  829. matrix_power(rshft, 2**100 + 2**10 + 2**5 + 2), self.rshft_2)
  830. assert_equal(
  831. matrix_power(rshft, 2**100 + 2**10 + 2**5 + 3), self.rshft_3)
  832. def test_power_is_zero(self, dt):
  833. def tz(M):
  834. mz = matrix_power(M, 0)
  835. assert_equal(mz, identity_like_generalized(M))
  836. assert_equal(mz.dtype, M.dtype)
  837. for mat in self.rshft_all:
  838. tz(mat.astype(dt))
  839. if dt != object:
  840. tz(self.stacked.astype(dt))
  841. def test_power_is_one(self, dt):
  842. def tz(mat):
  843. mz = matrix_power(mat, 1)
  844. assert_equal(mz, mat)
  845. assert_equal(mz.dtype, mat.dtype)
  846. for mat in self.rshft_all:
  847. tz(mat.astype(dt))
  848. if dt != object:
  849. tz(self.stacked.astype(dt))
  850. def test_power_is_two(self, dt):
  851. def tz(mat):
  852. mz = matrix_power(mat, 2)
  853. mmul = matmul if mat.dtype != object else dot
  854. assert_equal(mz, mmul(mat, mat))
  855. assert_equal(mz.dtype, mat.dtype)
  856. for mat in self.rshft_all:
  857. tz(mat.astype(dt))
  858. if dt != object:
  859. tz(self.stacked.astype(dt))
  860. def test_power_is_minus_one(self, dt):
  861. def tz(mat):
  862. invmat = matrix_power(mat, -1)
  863. mmul = matmul if mat.dtype != object else dot
  864. assert_almost_equal(
  865. mmul(invmat, mat), identity_like_generalized(mat))
  866. for mat in self.rshft_all:
  867. if dt not in self.dtnoinv:
  868. tz(mat.astype(dt))
  869. def test_exceptions_bad_power(self, dt):
  870. mat = self.rshft_0.astype(dt)
  871. assert_raises(TypeError, matrix_power, mat, 1.5)
  872. assert_raises(TypeError, matrix_power, mat, [1])
  873. def test_exceptions_non_square(self, dt):
  874. assert_raises(LinAlgError, matrix_power, np.array([1], dt), 1)
  875. assert_raises(LinAlgError, matrix_power, np.array([[1], [2]], dt), 1)
  876. assert_raises(LinAlgError, matrix_power, np.ones((4, 3, 2), dt), 1)
  877. def test_exceptions_not_invertible(self, dt):
  878. if dt in self.dtnoinv:
  879. return
  880. mat = self.noninv.astype(dt)
  881. assert_raises(LinAlgError, matrix_power, mat, -1)
  882. class TestEigvalshCases(HermitianTestCase, HermitianGeneralizedTestCase):
  883. def do(self, a, b, tags):
  884. # note that eigenvalue arrays returned by eig must be sorted since
  885. # their order isn't guaranteed.
  886. ev = linalg.eigvalsh(a, 'L')
  887. evalues, evectors = linalg.eig(a)
  888. evalues.sort(axis=-1)
  889. assert_allclose(ev, evalues, rtol=get_rtol(ev.dtype))
  890. ev2 = linalg.eigvalsh(a, 'U')
  891. assert_allclose(ev2, evalues, rtol=get_rtol(ev.dtype))
  892. class TestEigvalsh:
  893. @pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
  894. def test_types(self, dtype):
  895. x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
  896. w = np.linalg.eigvalsh(x)
  897. assert_equal(w.dtype, get_real_dtype(dtype))
  898. def test_invalid(self):
  899. x = np.array([[1, 0.5], [0.5, 1]], dtype=np.float32)
  900. assert_raises(ValueError, np.linalg.eigvalsh, x, UPLO="lrong")
  901. assert_raises(ValueError, np.linalg.eigvalsh, x, "lower")
  902. assert_raises(ValueError, np.linalg.eigvalsh, x, "upper")
  903. def test_UPLO(self):
  904. Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
  905. Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
  906. tgt = np.array([-1, 1], dtype=np.double)
  907. rtol = get_rtol(np.double)
  908. # Check default is 'L'
  909. w = np.linalg.eigvalsh(Klo)
  910. assert_allclose(w, tgt, rtol=rtol)
  911. # Check 'L'
  912. w = np.linalg.eigvalsh(Klo, UPLO='L')
  913. assert_allclose(w, tgt, rtol=rtol)
  914. # Check 'l'
  915. w = np.linalg.eigvalsh(Klo, UPLO='l')
  916. assert_allclose(w, tgt, rtol=rtol)
  917. # Check 'U'
  918. w = np.linalg.eigvalsh(Kup, UPLO='U')
  919. assert_allclose(w, tgt, rtol=rtol)
  920. # Check 'u'
  921. w = np.linalg.eigvalsh(Kup, UPLO='u')
  922. assert_allclose(w, tgt, rtol=rtol)
  923. def test_0_size(self):
  924. # Check that all kinds of 0-sized arrays work
  925. class ArraySubclass(np.ndarray):
  926. pass
  927. a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
  928. res = linalg.eigvalsh(a)
  929. assert_(res.dtype.type is np.float64)
  930. assert_equal((0, 1), res.shape)
  931. # This is just for documentation, it might make sense to change:
  932. assert_(isinstance(res, np.ndarray))
  933. a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
  934. res = linalg.eigvalsh(a)
  935. assert_(res.dtype.type is np.float32)
  936. assert_equal((0,), res.shape)
  937. # This is just for documentation, it might make sense to change:
  938. assert_(isinstance(res, np.ndarray))
  939. class TestEighCases(HermitianTestCase, HermitianGeneralizedTestCase):
  940. def do(self, a, b, tags):
  941. # note that eigenvalue arrays returned by eig must be sorted since
  942. # their order isn't guaranteed.
  943. ev, evc = linalg.eigh(a)
  944. evalues, evectors = linalg.eig(a)
  945. evalues.sort(axis=-1)
  946. assert_almost_equal(ev, evalues)
  947. assert_allclose(dot_generalized(a, evc),
  948. np.asarray(ev)[..., None, :] * np.asarray(evc),
  949. rtol=get_rtol(ev.dtype))
  950. ev2, evc2 = linalg.eigh(a, 'U')
  951. assert_almost_equal(ev2, evalues)
  952. assert_allclose(dot_generalized(a, evc2),
  953. np.asarray(ev2)[..., None, :] * np.asarray(evc2),
  954. rtol=get_rtol(ev.dtype), err_msg=repr(a))
  955. class TestEigh:
  956. @pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
  957. def test_types(self, dtype):
  958. x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
  959. w, v = np.linalg.eigh(x)
  960. assert_equal(w.dtype, get_real_dtype(dtype))
  961. assert_equal(v.dtype, dtype)
  962. def test_invalid(self):
  963. x = np.array([[1, 0.5], [0.5, 1]], dtype=np.float32)
  964. assert_raises(ValueError, np.linalg.eigh, x, UPLO="lrong")
  965. assert_raises(ValueError, np.linalg.eigh, x, "lower")
  966. assert_raises(ValueError, np.linalg.eigh, x, "upper")
  967. def test_UPLO(self):
  968. Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
  969. Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
  970. tgt = np.array([-1, 1], dtype=np.double)
  971. rtol = get_rtol(np.double)
  972. # Check default is 'L'
  973. w, v = np.linalg.eigh(Klo)
  974. assert_allclose(w, tgt, rtol=rtol)
  975. # Check 'L'
  976. w, v = np.linalg.eigh(Klo, UPLO='L')
  977. assert_allclose(w, tgt, rtol=rtol)
  978. # Check 'l'
  979. w, v = np.linalg.eigh(Klo, UPLO='l')
  980. assert_allclose(w, tgt, rtol=rtol)
  981. # Check 'U'
  982. w, v = np.linalg.eigh(Kup, UPLO='U')
  983. assert_allclose(w, tgt, rtol=rtol)
  984. # Check 'u'
  985. w, v = np.linalg.eigh(Kup, UPLO='u')
  986. assert_allclose(w, tgt, rtol=rtol)
  987. def test_0_size(self):
  988. # Check that all kinds of 0-sized arrays work
  989. class ArraySubclass(np.ndarray):
  990. pass
  991. a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
  992. res, res_v = linalg.eigh(a)
  993. assert_(res_v.dtype.type is np.float64)
  994. assert_(res.dtype.type is np.float64)
  995. assert_equal(a.shape, res_v.shape)
  996. assert_equal((0, 1), res.shape)
  997. # This is just for documentation, it might make sense to change:
  998. assert_(isinstance(a, np.ndarray))
  999. a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
  1000. res, res_v = linalg.eigh(a)
  1001. assert_(res_v.dtype.type is np.complex64)
  1002. assert_(res.dtype.type is np.float32)
  1003. assert_equal(a.shape, res_v.shape)
  1004. assert_equal((0,), res.shape)
  1005. # This is just for documentation, it might make sense to change:
  1006. assert_(isinstance(a, np.ndarray))
  1007. class _TestNormBase:
  1008. dt = None
  1009. dec = None
  1010. class _TestNormGeneral(_TestNormBase):
  1011. def test_empty(self):
  1012. assert_equal(norm([]), 0.0)
  1013. assert_equal(norm(array([], dtype=self.dt)), 0.0)
  1014. assert_equal(norm(atleast_2d(array([], dtype=self.dt))), 0.0)
  1015. def test_vector_return_type(self):
  1016. a = np.array([1, 0, 1])
  1017. exact_types = np.typecodes['AllInteger']
  1018. inexact_types = np.typecodes['AllFloat']
  1019. all_types = exact_types + inexact_types
  1020. for each_inexact_types in all_types:
  1021. at = a.astype(each_inexact_types)
  1022. an = norm(at, -np.inf)
  1023. assert_(issubclass(an.dtype.type, np.floating))
  1024. assert_almost_equal(an, 0.0)
  1025. with suppress_warnings() as sup:
  1026. sup.filter(RuntimeWarning, "divide by zero encountered")
  1027. an = norm(at, -1)
  1028. assert_(issubclass(an.dtype.type, np.floating))
  1029. assert_almost_equal(an, 0.0)
  1030. an = norm(at, 0)
  1031. assert_(issubclass(an.dtype.type, np.floating))
  1032. assert_almost_equal(an, 2)
  1033. an = norm(at, 1)
  1034. assert_(issubclass(an.dtype.type, np.floating))
  1035. assert_almost_equal(an, 2.0)
  1036. an = norm(at, 2)
  1037. assert_(issubclass(an.dtype.type, np.floating))
  1038. assert_almost_equal(an, an.dtype.type(2.0)**an.dtype.type(1.0/2.0))
  1039. an = norm(at, 4)
  1040. assert_(issubclass(an.dtype.type, np.floating))
  1041. assert_almost_equal(an, an.dtype.type(2.0)**an.dtype.type(1.0/4.0))
  1042. an = norm(at, np.inf)
  1043. assert_(issubclass(an.dtype.type, np.floating))
  1044. assert_almost_equal(an, 1.0)
  1045. def test_vector(self):
  1046. a = [1, 2, 3, 4]
  1047. b = [-1, -2, -3, -4]
  1048. c = [-1, 2, -3, 4]
  1049. def _test(v):
  1050. np.testing.assert_almost_equal(norm(v), 30 ** 0.5,
  1051. decimal=self.dec)
  1052. np.testing.assert_almost_equal(norm(v, inf), 4.0,
  1053. decimal=self.dec)
  1054. np.testing.assert_almost_equal(norm(v, -inf), 1.0,
  1055. decimal=self.dec)
  1056. np.testing.assert_almost_equal(norm(v, 1), 10.0,
  1057. decimal=self.dec)
  1058. np.testing.assert_almost_equal(norm(v, -1), 12.0 / 25,
  1059. decimal=self.dec)
  1060. np.testing.assert_almost_equal(norm(v, 2), 30 ** 0.5,
  1061. decimal=self.dec)
  1062. np.testing.assert_almost_equal(norm(v, -2), ((205. / 144) ** -0.5),
  1063. decimal=self.dec)
  1064. np.testing.assert_almost_equal(norm(v, 0), 4,
  1065. decimal=self.dec)
  1066. for v in (a, b, c,):
  1067. _test(v)
  1068. for v in (array(a, dtype=self.dt), array(b, dtype=self.dt),
  1069. array(c, dtype=self.dt)):
  1070. _test(v)
  1071. def test_axis(self):
  1072. # Vector norms.
  1073. # Compare the use of `axis` with computing the norm of each row
  1074. # or column separately.
  1075. A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt)
  1076. for order in [None, -1, 0, 1, 2, 3, np.Inf, -np.Inf]:
  1077. expected0 = [norm(A[:, k], ord=order) for k in range(A.shape[1])]
  1078. assert_almost_equal(norm(A, ord=order, axis=0), expected0)
  1079. expected1 = [norm(A[k, :], ord=order) for k in range(A.shape[0])]
  1080. assert_almost_equal(norm(A, ord=order, axis=1), expected1)
  1081. # Matrix norms.
  1082. B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)
  1083. nd = B.ndim
  1084. for order in [None, -2, 2, -1, 1, np.Inf, -np.Inf, 'fro']:
  1085. for axis in itertools.combinations(range(-nd, nd), 2):
  1086. row_axis, col_axis = axis
  1087. if row_axis < 0:
  1088. row_axis += nd
  1089. if col_axis < 0:
  1090. col_axis += nd
  1091. if row_axis == col_axis:
  1092. assert_raises(ValueError, norm, B, ord=order, axis=axis)
  1093. else:
  1094. n = norm(B, ord=order, axis=axis)
  1095. # The logic using k_index only works for nd = 3.
  1096. # This has to be changed if nd is increased.
  1097. k_index = nd - (row_axis + col_axis)
  1098. if row_axis < col_axis:
  1099. expected = [norm(B[:].take(k, axis=k_index), ord=order)
  1100. for k in range(B.shape[k_index])]
  1101. else:
  1102. expected = [norm(B[:].take(k, axis=k_index).T, ord=order)
  1103. for k in range(B.shape[k_index])]
  1104. assert_almost_equal(n, expected)
  1105. def test_keepdims(self):
  1106. A = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)
  1107. allclose_err = 'order {0}, axis = {1}'
  1108. shape_err = 'Shape mismatch found {0}, expected {1}, order={2}, axis={3}'
  1109. # check the order=None, axis=None case
  1110. expected = norm(A, ord=None, axis=None)
  1111. found = norm(A, ord=None, axis=None, keepdims=True)
  1112. assert_allclose(np.squeeze(found), expected,
  1113. err_msg=allclose_err.format(None, None))
  1114. expected_shape = (1, 1, 1)
  1115. assert_(found.shape == expected_shape,
  1116. shape_err.format(found.shape, expected_shape, None, None))
  1117. # Vector norms.
  1118. for order in [None, -1, 0, 1, 2, 3, np.Inf, -np.Inf]:
  1119. for k in range(A.ndim):
  1120. expected = norm(A, ord=order, axis=k)
  1121. found = norm(A, ord=order, axis=k, keepdims=True)
  1122. assert_allclose(np.squeeze(found), expected,
  1123. err_msg=allclose_err.format(order, k))
  1124. expected_shape = list(A.shape)
  1125. expected_shape[k] = 1
  1126. expected_shape = tuple(expected_shape)
  1127. assert_(found.shape == expected_shape,
  1128. shape_err.format(found.shape, expected_shape, order, k))
  1129. # Matrix norms.
  1130. for order in [None, -2, 2, -1, 1, np.Inf, -np.Inf, 'fro', 'nuc']:
  1131. for k in itertools.permutations(range(A.ndim), 2):
  1132. expected = norm(A, ord=order, axis=k)
  1133. found = norm(A, ord=order, axis=k, keepdims=True)
  1134. assert_allclose(np.squeeze(found), expected,
  1135. err_msg=allclose_err.format(order, k))
  1136. expected_shape = list(A.shape)
  1137. expected_shape[k[0]] = 1
  1138. expected_shape[k[1]] = 1
  1139. expected_shape = tuple(expected_shape)
  1140. assert_(found.shape == expected_shape,
  1141. shape_err.format(found.shape, expected_shape, order, k))
  1142. class _TestNorm2D(_TestNormBase):
  1143. # Define the part for 2d arrays separately, so we can subclass this
  1144. # and run the tests using np.matrix in matrixlib.tests.test_matrix_linalg.
  1145. array = np.array
  1146. def test_matrix_empty(self):
  1147. assert_equal(norm(self.array([[]], dtype=self.dt)), 0.0)
  1148. def test_matrix_return_type(self):
  1149. a = self.array([[1, 0, 1], [0, 1, 1]])
  1150. exact_types = np.typecodes['AllInteger']
  1151. # float32, complex64, float64, complex128 types are the only types
  1152. # allowed by `linalg`, which performs the matrix operations used
  1153. # within `norm`.
  1154. inexact_types = 'fdFD'
  1155. all_types = exact_types + inexact_types
  1156. for each_inexact_types in all_types:
  1157. at = a.astype(each_inexact_types)
  1158. an = norm(at, -np.inf)
  1159. assert_(issubclass(an.dtype.type, np.floating))
  1160. assert_almost_equal(an, 2.0)
  1161. with suppress_warnings() as sup:
  1162. sup.filter(RuntimeWarning, "divide by zero encountered")
  1163. an = norm(at, -1)
  1164. assert_(issubclass(an.dtype.type, np.floating))
  1165. assert_almost_equal(an, 1.0)
  1166. an = norm(at, 1)
  1167. assert_(issubclass(an.dtype.type, np.floating))
  1168. assert_almost_equal(an, 2.0)
  1169. an = norm(at, 2)
  1170. assert_(issubclass(an.dtype.type, np.floating))
  1171. assert_almost_equal(an, 3.0**(1.0/2.0))
  1172. an = norm(at, -2)
  1173. assert_(issubclass(an.dtype.type, np.floating))
  1174. assert_almost_equal(an, 1.0)
  1175. an = norm(at, np.inf)
  1176. assert_(issubclass(an.dtype.type, np.floating))
  1177. assert_almost_equal(an, 2.0)
  1178. an = norm(at, 'fro')
  1179. assert_(issubclass(an.dtype.type, np.floating))
  1180. assert_almost_equal(an, 2.0)
  1181. an = norm(at, 'nuc')
  1182. assert_(issubclass(an.dtype.type, np.floating))
  1183. # Lower bar needed to support low precision floats.
  1184. # They end up being off by 1 in the 7th place.
  1185. np.testing.assert_almost_equal(an, 2.7320508075688772, decimal=6)
  1186. def test_matrix_2x2(self):
  1187. A = self.array([[1, 3], [5, 7]], dtype=self.dt)
  1188. assert_almost_equal(norm(A), 84 ** 0.5)
  1189. assert_almost_equal(norm(A, 'fro'), 84 ** 0.5)
  1190. assert_almost_equal(norm(A, 'nuc'), 10.0)
  1191. assert_almost_equal(norm(A, inf), 12.0)
  1192. assert_almost_equal(norm(A, -inf), 4.0)
  1193. assert_almost_equal(norm(A, 1), 10.0)
  1194. assert_almost_equal(norm(A, -1), 6.0)
  1195. assert_almost_equal(norm(A, 2), 9.1231056256176615)
  1196. assert_almost_equal(norm(A, -2), 0.87689437438234041)
  1197. assert_raises(ValueError, norm, A, 'nofro')
  1198. assert_raises(ValueError, norm, A, -3)
  1199. assert_raises(ValueError, norm, A, 0)
  1200. def test_matrix_3x3(self):
  1201. # This test has been added because the 2x2 example
  1202. # happened to have equal nuclear norm and induced 1-norm.
  1203. # The 1/10 scaling factor accommodates the absolute tolerance
  1204. # used in assert_almost_equal.
  1205. A = (1 / 10) * \
  1206. self.array([[1, 2, 3], [6, 0, 5], [3, 2, 1]], dtype=self.dt)
  1207. assert_almost_equal(norm(A), (1 / 10) * 89 ** 0.5)
  1208. assert_almost_equal(norm(A, 'fro'), (1 / 10) * 89 ** 0.5)
  1209. assert_almost_equal(norm(A, 'nuc'), 1.3366836911774836)
  1210. assert_almost_equal(norm(A, inf), 1.1)
  1211. assert_almost_equal(norm(A, -inf), 0.6)
  1212. assert_almost_equal(norm(A, 1), 1.0)
  1213. assert_almost_equal(norm(A, -1), 0.4)
  1214. assert_almost_equal(norm(A, 2), 0.88722940323461277)
  1215. assert_almost_equal(norm(A, -2), 0.19456584790481812)
  1216. def test_bad_args(self):
  1217. # Check that bad arguments raise the appropriate exceptions.
  1218. A = self.array([[1, 2, 3], [4, 5, 6]], dtype=self.dt)
  1219. B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)
  1220. # Using `axis=<integer>` or passing in a 1-D array implies vector
  1221. # norms are being computed, so also using `ord='fro'`
  1222. # or `ord='nuc'` or any other string raises a ValueError.
  1223. assert_raises(ValueError, norm, A, 'fro', 0)
  1224. assert_raises(ValueError, norm, A, 'nuc', 0)
  1225. assert_raises(ValueError, norm, [3, 4], 'fro', None)
  1226. assert_raises(ValueError, norm, [3, 4], 'nuc', None)
  1227. assert_raises(ValueError, norm, [3, 4], 'test', None)
  1228. # Similarly, norm should raise an exception when ord is any finite
  1229. # number other than 1, 2, -1 or -2 when computing matrix norms.
  1230. for order in [0, 3]:
  1231. assert_raises(ValueError, norm, A, order, None)
  1232. assert_raises(ValueError, norm, A, order, (0, 1))
  1233. assert_raises(ValueError, norm, B, order, (1, 2))
  1234. # Invalid axis
  1235. assert_raises(np.AxisError, norm, B, None, 3)
  1236. assert_raises(np.AxisError, norm, B, None, (2, 3))
  1237. assert_raises(ValueError, norm, B, None, (0, 1, 2))
  1238. class _TestNorm(_TestNorm2D, _TestNormGeneral):
  1239. pass
  1240. class TestNorm_NonSystematic:
  1241. def test_longdouble_norm(self):
  1242. # Non-regression test: p-norm of longdouble would previously raise
  1243. # UnboundLocalError.
  1244. x = np.arange(10, dtype=np.longdouble)
  1245. old_assert_almost_equal(norm(x, ord=3), 12.65, decimal=2)
  1246. def test_intmin(self):
  1247. # Non-regression test: p-norm of signed integer would previously do
  1248. # float cast and abs in the wrong order.
  1249. x = np.array([-2 ** 31], dtype=np.int32)
  1250. old_assert_almost_equal(norm(x, ord=3), 2 ** 31, decimal=5)
  1251. def test_complex_high_ord(self):
  1252. # gh-4156
  1253. d = np.empty((2,), dtype=np.clongdouble)
  1254. d[0] = 6 + 7j
  1255. d[1] = -6 + 7j
  1256. res = 11.615898132184
  1257. old_assert_almost_equal(np.linalg.norm(d, ord=3), res, decimal=10)
  1258. d = d.astype(np.complex128)
  1259. old_assert_almost_equal(np.linalg.norm(d, ord=3), res, decimal=9)
  1260. d = d.astype(np.complex64)
  1261. old_assert_almost_equal(np.linalg.norm(d, ord=3), res, decimal=5)
  1262. # Separate definitions so we can use them for matrix tests.
  1263. class _TestNormDoubleBase(_TestNormBase):
  1264. dt = np.double
  1265. dec = 12
  1266. class _TestNormSingleBase(_TestNormBase):
  1267. dt = np.float32
  1268. dec = 6
  1269. class _TestNormInt64Base(_TestNormBase):
  1270. dt = np.int64
  1271. dec = 12
  1272. class TestNormDouble(_TestNorm, _TestNormDoubleBase):
  1273. pass
  1274. class TestNormSingle(_TestNorm, _TestNormSingleBase):
  1275. pass
  1276. class TestNormInt64(_TestNorm, _TestNormInt64Base):
  1277. pass
  1278. class TestMatrixRank:
  1279. def test_matrix_rank(self):
  1280. # Full rank matrix
  1281. assert_equal(4, matrix_rank(np.eye(4)))
  1282. # rank deficient matrix
  1283. I = np.eye(4)
  1284. I[-1, -1] = 0.
  1285. assert_equal(matrix_rank(I), 3)
  1286. # All zeros - zero rank
  1287. assert_equal(matrix_rank(np.zeros((4, 4))), 0)
  1288. # 1 dimension - rank 1 unless all 0
  1289. assert_equal(matrix_rank([1, 0, 0, 0]), 1)
  1290. assert_equal(matrix_rank(np.zeros((4,))), 0)
  1291. # accepts array-like
  1292. assert_equal(matrix_rank([1]), 1)
  1293. # greater than 2 dimensions treated as stacked matrices
  1294. ms = np.array([I, np.eye(4), np.zeros((4,4))])
  1295. assert_equal(matrix_rank(ms), np.array([3, 4, 0]))
  1296. # works on scalar
  1297. assert_equal(matrix_rank(1), 1)
  1298. def test_symmetric_rank(self):
  1299. assert_equal(4, matrix_rank(np.eye(4), hermitian=True))
  1300. assert_equal(1, matrix_rank(np.ones((4, 4)), hermitian=True))
  1301. assert_equal(0, matrix_rank(np.zeros((4, 4)), hermitian=True))
  1302. # rank deficient matrix
  1303. I = np.eye(4)
  1304. I[-1, -1] = 0.
  1305. assert_equal(3, matrix_rank(I, hermitian=True))
  1306. # manually supplied tolerance
  1307. I[-1, -1] = 1e-8
  1308. assert_equal(4, matrix_rank(I, hermitian=True, tol=0.99e-8))
  1309. assert_equal(3, matrix_rank(I, hermitian=True, tol=1.01e-8))
  1310. def test_reduced_rank():
  1311. # Test matrices with reduced rank
  1312. rng = np.random.RandomState(20120714)
  1313. for i in range(100):
  1314. # Make a rank deficient matrix
  1315. X = rng.normal(size=(40, 10))
  1316. X[:, 0] = X[:, 1] + X[:, 2]
  1317. # Assert that matrix_rank detected deficiency
  1318. assert_equal(matrix_rank(X), 9)
  1319. X[:, 3] = X[:, 4] + X[:, 5]
  1320. assert_equal(matrix_rank(X), 8)
  1321. class TestQR:
  1322. # Define the array class here, so run this on matrices elsewhere.
  1323. array = np.array
  1324. def check_qr(self, a):
  1325. # This test expects the argument `a` to be an ndarray or
  1326. # a subclass of an ndarray of inexact type.
  1327. a_type = type(a)
  1328. a_dtype = a.dtype
  1329. m, n = a.shape
  1330. k = min(m, n)
  1331. # mode == 'complete'
  1332. q, r = linalg.qr(a, mode='complete')
  1333. assert_(q.dtype == a_dtype)
  1334. assert_(r.dtype == a_dtype)
  1335. assert_(isinstance(q, a_type))
  1336. assert_(isinstance(r, a_type))
  1337. assert_(q.shape == (m, m))
  1338. assert_(r.shape == (m, n))
  1339. assert_almost_equal(dot(q, r), a)
  1340. assert_almost_equal(dot(q.T.conj(), q), np.eye(m))
  1341. assert_almost_equal(np.triu(r), r)
  1342. # mode == 'reduced'
  1343. q1, r1 = linalg.qr(a, mode='reduced')
  1344. assert_(q1.dtype == a_dtype)
  1345. assert_(r1.dtype == a_dtype)
  1346. assert_(isinstance(q1, a_type))
  1347. assert_(isinstance(r1, a_type))
  1348. assert_(q1.shape == (m, k))
  1349. assert_(r1.shape == (k, n))
  1350. assert_almost_equal(dot(q1, r1), a)
  1351. assert_almost_equal(dot(q1.T.conj(), q1), np.eye(k))
  1352. assert_almost_equal(np.triu(r1), r1)
  1353. # mode == 'r'
  1354. r2 = linalg.qr(a, mode='r')
  1355. assert_(r2.dtype == a_dtype)
  1356. assert_(isinstance(r2, a_type))
  1357. assert_almost_equal(r2, r1)
  1358. @pytest.mark.parametrize(["m", "n"], [
  1359. (3, 0),
  1360. (0, 3),
  1361. (0, 0)
  1362. ])
  1363. def test_qr_empty(self, m, n):
  1364. k = min(m, n)
  1365. a = np.empty((m, n))
  1366. self.check_qr(a)
  1367. h, tau = np.linalg.qr(a, mode='raw')
  1368. assert_equal(h.dtype, np.double)
  1369. assert_equal(tau.dtype, np.double)
  1370. assert_equal(h.shape, (n, m))
  1371. assert_equal(tau.shape, (k,))
  1372. def test_mode_raw(self):
  1373. # The factorization is not unique and varies between libraries,
  1374. # so it is not possible to check against known values. Functional
  1375. # testing is a possibility, but awaits the exposure of more
  1376. # of the functions in lapack_lite. Consequently, this test is
  1377. # very limited in scope. Note that the results are in FORTRAN
  1378. # order, hence the h arrays are transposed.
  1379. a = self.array([[1, 2], [3, 4], [5, 6]], dtype=np.double)
  1380. # Test double
  1381. h, tau = linalg.qr(a, mode='raw')
  1382. assert_(h.dtype == np.double)
  1383. assert_(tau.dtype == np.double)
  1384. assert_(h.shape == (2, 3))
  1385. assert_(tau.shape == (2,))
  1386. h, tau = linalg.qr(a.T, mode='raw')
  1387. assert_(h.dtype == np.double)
  1388. assert_(tau.dtype == np.double)
  1389. assert_(h.shape == (3, 2))
  1390. assert_(tau.shape == (2,))
  1391. def test_mode_all_but_economic(self):
  1392. a = self.array([[1, 2], [3, 4]])
  1393. b = self.array([[1, 2], [3, 4], [5, 6]])
  1394. for dt in "fd":
  1395. m1 = a.astype(dt)
  1396. m2 = b.astype(dt)
  1397. self.check_qr(m1)
  1398. self.check_qr(m2)
  1399. self.check_qr(m2.T)
  1400. for dt in "fd":
  1401. m1 = 1 + 1j * a.astype(dt)
  1402. m2 = 1 + 1j * b.astype(dt)
  1403. self.check_qr(m1)
  1404. self.check_qr(m2)
  1405. self.check_qr(m2.T)
  1406. class TestCholesky:
  1407. # TODO: are there no other tests for cholesky?
  1408. def test_basic_property(self):
  1409. # Check A = L L^H
  1410. shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
  1411. dtypes = (np.float32, np.float64, np.complex64, np.complex128)
  1412. for shape, dtype in itertools.product(shapes, dtypes):
  1413. np.random.seed(1)
  1414. a = np.random.randn(*shape)
  1415. if np.issubdtype(dtype, np.complexfloating):
  1416. a = a + 1j*np.random.randn(*shape)
  1417. t = list(range(len(shape)))
  1418. t[-2:] = -1, -2
  1419. a = np.matmul(a.transpose(t).conj(), a)
  1420. a = np.asarray(a, dtype=dtype)
  1421. c = np.linalg.cholesky(a)
  1422. b = np.matmul(c, c.transpose(t).conj())
  1423. assert_allclose(b, a,
  1424. err_msg=f'{shape} {dtype}\n{a}\n{c}',
  1425. atol=500 * a.shape[0] * np.finfo(dtype).eps)
  1426. def test_0_size(self):
  1427. class ArraySubclass(np.ndarray):
  1428. pass
  1429. a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
  1430. res = linalg.cholesky(a)
  1431. assert_equal(a.shape, res.shape)
  1432. assert_(res.dtype.type is np.float64)
  1433. # for documentation purpose:
  1434. assert_(isinstance(res, np.ndarray))
  1435. a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
  1436. res = linalg.cholesky(a)
  1437. assert_equal(a.shape, res.shape)
  1438. assert_(res.dtype.type is np.complex64)
  1439. assert_(isinstance(res, np.ndarray))
  1440. def test_byteorder_check():
  1441. # Byte order check should pass for native order
  1442. if sys.byteorder == 'little':
  1443. native = '<'
  1444. else:
  1445. native = '>'
  1446. for dtt in (np.float32, np.float64):
  1447. arr = np.eye(4, dtype=dtt)
  1448. n_arr = arr.newbyteorder(native)
  1449. sw_arr = arr.newbyteorder('S').byteswap()
  1450. assert_equal(arr.dtype.byteorder, '=')
  1451. for routine in (linalg.inv, linalg.det, linalg.pinv):
  1452. # Normal call
  1453. res = routine(arr)
  1454. # Native but not '='
  1455. assert_array_equal(res, routine(n_arr))
  1456. # Swapped
  1457. assert_array_equal(res, routine(sw_arr))
  1458. def test_generalized_raise_multiloop():
  1459. # It should raise an error even if the error doesn't occur in the
  1460. # last iteration of the ufunc inner loop
  1461. invertible = np.array([[1, 2], [3, 4]])
  1462. non_invertible = np.array([[1, 1], [1, 1]])
  1463. x = np.zeros([4, 4, 2, 2])[1::2]
  1464. x[...] = invertible
  1465. x[0, 0] = non_invertible
  1466. assert_raises(np.linalg.LinAlgError, np.linalg.inv, x)
  1467. def test_xerbla_override():
  1468. # Check that our xerbla has been successfully linked in. If it is not,
  1469. # the default xerbla routine is called, which prints a message to stdout
  1470. # and may, or may not, abort the process depending on the LAPACK package.
  1471. XERBLA_OK = 255
  1472. try:
  1473. pid = os.fork()
  1474. except (OSError, AttributeError):
  1475. # fork failed, or not running on POSIX
  1476. pytest.skip("Not POSIX or fork failed.")
  1477. if pid == 0:
  1478. # child; close i/o file handles
  1479. os.close(1)
  1480. os.close(0)
  1481. # Avoid producing core files.
  1482. import resource
  1483. resource.setrlimit(resource.RLIMIT_CORE, (0, 0))
  1484. # These calls may abort.
  1485. try:
  1486. np.linalg.lapack_lite.xerbla()
  1487. except ValueError:
  1488. pass
  1489. except Exception:
  1490. os._exit(os.EX_CONFIG)
  1491. try:
  1492. a = np.array([[1.]])
  1493. np.linalg.lapack_lite.dorgqr(
  1494. 1, 1, 1, a,
  1495. 0, # <- invalid value
  1496. a, a, 0, 0)
  1497. except ValueError as e:
  1498. if "DORGQR parameter number 5" in str(e):
  1499. # success, reuse error code to mark success as
  1500. # FORTRAN STOP returns as success.
  1501. os._exit(XERBLA_OK)
  1502. # Did not abort, but our xerbla was not linked in.
  1503. os._exit(os.EX_CONFIG)
  1504. else:
  1505. # parent
  1506. pid, status = os.wait()
  1507. if os.WEXITSTATUS(status) != XERBLA_OK:
  1508. pytest.skip('Numpy xerbla not linked in.')
  1509. @pytest.mark.slow
  1510. def test_sdot_bug_8577():
  1511. # Regression test that loading certain other libraries does not
  1512. # result to wrong results in float32 linear algebra.
  1513. #
  1514. # There's a bug gh-8577 on OSX that can trigger this, and perhaps
  1515. # there are also other situations in which it occurs.
  1516. #
  1517. # Do the check in a separate process.
  1518. bad_libs = ['PyQt5.QtWidgets', 'IPython']
  1519. template = textwrap.dedent("""
  1520. import sys
  1521. {before}
  1522. try:
  1523. import {bad_lib}
  1524. except ImportError:
  1525. sys.exit(0)
  1526. {after}
  1527. x = np.ones(2, dtype=np.float32)
  1528. sys.exit(0 if np.allclose(x.dot(x), 2.0) else 1)
  1529. """)
  1530. for bad_lib in bad_libs:
  1531. code = template.format(before="import numpy as np", after="",
  1532. bad_lib=bad_lib)
  1533. subprocess.check_call([sys.executable, "-c", code])
  1534. # Swapped import order
  1535. code = template.format(after="import numpy as np", before="",
  1536. bad_lib=bad_lib)
  1537. subprocess.check_call([sys.executable, "-c", code])
  1538. class TestMultiDot:
  1539. def test_basic_function_with_three_arguments(self):
  1540. # multi_dot with three arguments uses a fast hand coded algorithm to
  1541. # determine the optimal order. Therefore test it separately.
  1542. A = np.random.random((6, 2))
  1543. B = np.random.random((2, 6))
  1544. C = np.random.random((6, 2))
  1545. assert_almost_equal(multi_dot([A, B, C]), A.dot(B).dot(C))
  1546. assert_almost_equal(multi_dot([A, B, C]), np.dot(A, np.dot(B, C)))
  1547. def test_basic_function_with_two_arguments(self):
  1548. # separate code path with two arguments
  1549. A = np.random.random((6, 2))
  1550. B = np.random.random((2, 6))
  1551. assert_almost_equal(multi_dot([A, B]), A.dot(B))
  1552. assert_almost_equal(multi_dot([A, B]), np.dot(A, B))
  1553. def test_basic_function_with_dynamic_programing_optimization(self):
  1554. # multi_dot with four or more arguments uses the dynamic programing
  1555. # optimization and therefore deserve a separate
  1556. A = np.random.random((6, 2))
  1557. B = np.random.random((2, 6))
  1558. C = np.random.random((6, 2))
  1559. D = np.random.random((2, 1))
  1560. assert_almost_equal(multi_dot([A, B, C, D]), A.dot(B).dot(C).dot(D))
  1561. def test_vector_as_first_argument(self):
  1562. # The first argument can be 1-D
  1563. A1d = np.random.random(2) # 1-D
  1564. B = np.random.random((2, 6))
  1565. C = np.random.random((6, 2))
  1566. D = np.random.random((2, 2))
  1567. # the result should be 1-D
  1568. assert_equal(multi_dot([A1d, B, C, D]).shape, (2,))
  1569. def test_vector_as_last_argument(self):
  1570. # The last argument can be 1-D
  1571. A = np.random.random((6, 2))
  1572. B = np.random.random((2, 6))
  1573. C = np.random.random((6, 2))
  1574. D1d = np.random.random(2) # 1-D
  1575. # the result should be 1-D
  1576. assert_equal(multi_dot([A, B, C, D1d]).shape, (6,))
  1577. def test_vector_as_first_and_last_argument(self):
  1578. # The first and last arguments can be 1-D
  1579. A1d = np.random.random(2) # 1-D
  1580. B = np.random.random((2, 6))
  1581. C = np.random.random((6, 2))
  1582. D1d = np.random.random(2) # 1-D
  1583. # the result should be a scalar
  1584. assert_equal(multi_dot([A1d, B, C, D1d]).shape, ())
  1585. def test_three_arguments_and_out(self):
  1586. # multi_dot with three arguments uses a fast hand coded algorithm to
  1587. # determine the optimal order. Therefore test it separately.
  1588. A = np.random.random((6, 2))
  1589. B = np.random.random((2, 6))
  1590. C = np.random.random((6, 2))
  1591. out = np.zeros((6, 2))
  1592. ret = multi_dot([A, B, C], out=out)
  1593. assert out is ret
  1594. assert_almost_equal(out, A.dot(B).dot(C))
  1595. assert_almost_equal(out, np.dot(A, np.dot(B, C)))
  1596. def test_two_arguments_and_out(self):
  1597. # separate code path with two arguments
  1598. A = np.random.random((6, 2))
  1599. B = np.random.random((2, 6))
  1600. out = np.zeros((6, 6))
  1601. ret = multi_dot([A, B], out=out)
  1602. assert out is ret
  1603. assert_almost_equal(out, A.dot(B))
  1604. assert_almost_equal(out, np.dot(A, B))
  1605. def test_dynamic_programing_optimization_and_out(self):
  1606. # multi_dot with four or more arguments uses the dynamic programing
  1607. # optimization and therefore deserve a separate test
  1608. A = np.random.random((6, 2))
  1609. B = np.random.random((2, 6))
  1610. C = np.random.random((6, 2))
  1611. D = np.random.random((2, 1))
  1612. out = np.zeros((6, 1))
  1613. ret = multi_dot([A, B, C, D], out=out)
  1614. assert out is ret
  1615. assert_almost_equal(out, A.dot(B).dot(C).dot(D))
  1616. def test_dynamic_programming_logic(self):
  1617. # Test for the dynamic programming part
  1618. # This test is directly taken from Cormen page 376.
  1619. arrays = [np.random.random((30, 35)),
  1620. np.random.random((35, 15)),
  1621. np.random.random((15, 5)),
  1622. np.random.random((5, 10)),
  1623. np.random.random((10, 20)),
  1624. np.random.random((20, 25))]
  1625. m_expected = np.array([[0., 15750., 7875., 9375., 11875., 15125.],
  1626. [0., 0., 2625., 4375., 7125., 10500.],
  1627. [0., 0., 0., 750., 2500., 5375.],
  1628. [0., 0., 0., 0., 1000., 3500.],
  1629. [0., 0., 0., 0., 0., 5000.],
  1630. [0., 0., 0., 0., 0., 0.]])
  1631. s_expected = np.array([[0, 1, 1, 3, 3, 3],
  1632. [0, 0, 2, 3, 3, 3],
  1633. [0, 0, 0, 3, 3, 3],
  1634. [0, 0, 0, 0, 4, 5],
  1635. [0, 0, 0, 0, 0, 5],
  1636. [0, 0, 0, 0, 0, 0]], dtype=int)
  1637. s_expected -= 1 # Cormen uses 1-based index, python does not.
  1638. s, m = _multi_dot_matrix_chain_order(arrays, return_costs=True)
  1639. # Only the upper triangular part (without the diagonal) is interesting.
  1640. assert_almost_equal(np.triu(s[:-1, 1:]),
  1641. np.triu(s_expected[:-1, 1:]))
  1642. assert_almost_equal(np.triu(m), np.triu(m_expected))
  1643. def test_too_few_input_arrays(self):
  1644. assert_raises(ValueError, multi_dot, [])
  1645. assert_raises(ValueError, multi_dot, [np.random.random((3, 3))])
  1646. class TestTensorinv:
  1647. @pytest.mark.parametrize("arr, ind", [
  1648. (np.ones((4, 6, 8, 2)), 2),
  1649. (np.ones((3, 3, 2)), 1),
  1650. ])
  1651. def test_non_square_handling(self, arr, ind):
  1652. with assert_raises(LinAlgError):
  1653. linalg.tensorinv(arr, ind=ind)
  1654. @pytest.mark.parametrize("shape, ind", [
  1655. # examples from docstring
  1656. ((4, 6, 8, 3), 2),
  1657. ((24, 8, 3), 1),
  1658. ])
  1659. def test_tensorinv_shape(self, shape, ind):
  1660. a = np.eye(24)
  1661. a.shape = shape
  1662. ainv = linalg.tensorinv(a=a, ind=ind)
  1663. expected = a.shape[ind:] + a.shape[:ind]
  1664. actual = ainv.shape
  1665. assert_equal(actual, expected)
  1666. @pytest.mark.parametrize("ind", [
  1667. 0, -2,
  1668. ])
  1669. def test_tensorinv_ind_limit(self, ind):
  1670. a = np.eye(24)
  1671. a.shape = (4, 6, 8, 3)
  1672. with assert_raises(ValueError):
  1673. linalg.tensorinv(a=a, ind=ind)
  1674. def test_tensorinv_result(self):
  1675. # mimic a docstring example
  1676. a = np.eye(24)
  1677. a.shape = (24, 8, 3)
  1678. ainv = linalg.tensorinv(a, ind=1)
  1679. b = np.ones(24)
  1680. assert_allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
  1681. def test_unsupported_commontype():
  1682. # linalg gracefully handles unsupported type
  1683. arr = np.array([[1, -2], [2, 5]], dtype='float16')
  1684. with assert_raises_regex(TypeError, "unsupported in linalg"):
  1685. linalg.cholesky(arr)
  1686. @pytest.mark.slow
  1687. @pytest.mark.xfail(not HAS_LAPACK64, run=False,
  1688. reason="Numpy not compiled with 64-bit BLAS/LAPACK")
  1689. @requires_memory(free_bytes=16e9)
  1690. def test_blas64_dot():
  1691. n = 2**32
  1692. a = np.zeros([1, n], dtype=np.float32)
  1693. b = np.ones([1, 1], dtype=np.float32)
  1694. a[0,-1] = 1
  1695. c = np.dot(b, a)
  1696. assert_equal(c[0,-1], 1)
  1697. @pytest.mark.xfail(not HAS_LAPACK64,
  1698. reason="Numpy not compiled with 64-bit BLAS/LAPACK")
  1699. def test_blas64_geqrf_lwork_smoketest():
  1700. # Smoke test LAPACK geqrf lwork call with 64-bit integers
  1701. dtype = np.float64
  1702. lapack_routine = np.linalg.lapack_lite.dgeqrf
  1703. m = 2**32 + 1
  1704. n = 2**32 + 1
  1705. lda = m
  1706. # Dummy arrays, not referenced by the lapack routine, so don't
  1707. # need to be of the right size
  1708. a = np.zeros([1, 1], dtype=dtype)
  1709. work = np.zeros([1], dtype=dtype)
  1710. tau = np.zeros([1], dtype=dtype)
  1711. # Size query
  1712. results = lapack_routine(m, n, a, lda, tau, work, -1, 0)
  1713. assert_equal(results['info'], 0)
  1714. assert_equal(results['m'], m)
  1715. assert_equal(results['n'], m)
  1716. # Should result to an integer of a reasonable size
  1717. lwork = int(work.item())
  1718. assert_(2**32 < lwork < 2**42)