import itertools
import platform
import sys
import numpy as np
from numpy.testing import (assert_equal, assert_almost_equal,
assert_array_almost_equal, assert_array_equal,
assert_, assert_allclose)
import pytest
from pytest import raises as assert_raises
from scipy.linalg import (eig, eigvals, lu, svd, svdvals, cholesky, qr,
schur, rsf2csf, lu_solve, lu_factor, solve, diagsvd,
hessenberg, rq, eig_banded, eigvals_banded, eigh,
eigvalsh, qr_multiply, qz, orth, ordqz,
subspace_angles, hadamard, eigvalsh_tridiagonal,
eigh_tridiagonal, null_space, cdf2rdf, LinAlgError)
from scipy.linalg.lapack import (dgbtrf, dgbtrs, zgbtrf, zgbtrs, dsbev,
dsbevd, dsbevx, zhbevd, zhbevx)
from scipy.linalg._misc import norm
from scipy.linalg._decomp_qz import _select_function
from scipy.stats import ortho_group
from numpy import (array, diag, full, linalg, argsort, zeros, arange,
float32, complex64, ravel, sqrt, iscomplex, shape, sort,
sign, asarray, isfinite, ndarray, eye,)
from scipy.linalg._testutils import assert_no_overwrite
from scipy.sparse._sputils import matrix
from scipy._lib._testutils import check_free_memory
from scipy.linalg.blas import HAS_ILP64
try:
from scipy.__config__ import CONFIG
except ImportError:
CONFIG = None
def _random_hermitian_matrix(n, posdef=False, dtype=float):
"Generate random sym/hermitian array of the given size n"
if dtype in COMPLEX_DTYPES:
A = np.random.rand(n, n) + np.random.rand(n, n)*1.0j
A = (A + A.conj().T)/2
else:
A = np.random.rand(n, n)
A = (A + A.T)/2
if posdef:
A += sqrt(2*n)*np.eye(n)
return A.astype(dtype)
REAL_DTYPES = [np.float32, np.float64]
COMPLEX_DTYPES = [np.complex64, np.complex128]
DTYPES = REAL_DTYPES + COMPLEX_DTYPES
# XXX: This function should not be defined here, but somewhere in
# scipy.linalg namespace
def symrand(dim_or_eigv, rng):
"""Return a random symmetric (Hermitian) matrix.
If 'dim_or_eigv' is an integer N, return a NxN matrix, with eigenvalues
uniformly distributed on (-1,1).
If 'dim_or_eigv' is 1-D real array 'a', return a matrix whose
eigenvalues are 'a'.
"""
if isinstance(dim_or_eigv, int):
dim = dim_or_eigv
d = rng.random(dim)*2 - 1
elif (isinstance(dim_or_eigv, ndarray) and
len(dim_or_eigv.shape) == 1):
dim = dim_or_eigv.shape[0]
d = dim_or_eigv
else:
raise TypeError("input type not supported.")
v = ortho_group.rvs(dim)
h = v.T.conj() @ diag(d) @ v
# to avoid roundoff errors, symmetrize the matrix (again)
h = 0.5*(h.T+h)
return h
class TestEigVals:
def test_simple(self):
a = [[1, 2, 3], [1, 2, 3], [2, 5, 6]]
w = eigvals(a)
exact_w = [(9+sqrt(93))/2, 0, (9-sqrt(93))/2]
assert_array_almost_equal(w, exact_w)
def test_simple_tr(self):
a = array([[1, 2, 3], [1, 2, 3], [2, 5, 6]], 'd').T
a = a.copy()
a = a.T
w = eigvals(a)
exact_w = [(9+sqrt(93))/2, 0, (9-sqrt(93))/2]
assert_array_almost_equal(w, exact_w)
def test_simple_complex(self):
a = [[1, 2, 3], [1, 2, 3], [2, 5, 6+1j]]
w = eigvals(a)
exact_w = [(9+1j+sqrt(92+6j))/2,
0,
(9+1j-sqrt(92+6j))/2]
assert_array_almost_equal(w, exact_w)
def test_finite(self):
a = [[1, 2, 3], [1, 2, 3], [2, 5, 6]]
w = eigvals(a, check_finite=False)
exact_w = [(9+sqrt(93))/2, 0, (9-sqrt(93))/2]
assert_array_almost_equal(w, exact_w)
class TestEig:
def test_simple(self):
a = array([[1, 2, 3], [1, 2, 3], [2, 5, 6]])
w, v = eig(a)
exact_w = [(9+sqrt(93))/2, 0, (9-sqrt(93))/2]
v0 = array([1, 1, (1+sqrt(93)/3)/2])
v1 = array([3., 0, -1])
v2 = array([1, 1, (1-sqrt(93)/3)/2])
v0 = v0 / norm(v0)
v1 = v1 / norm(v1)
v2 = v2 / norm(v2)
assert_array_almost_equal(w, exact_w)
assert_array_almost_equal(v0, v[:, 0]*sign(v[0, 0]))
assert_array_almost_equal(v1, v[:, 1]*sign(v[0, 1]))
assert_array_almost_equal(v2, v[:, 2]*sign(v[0, 2]))
for i in range(3):
assert_array_almost_equal(a @ v[:, i], w[i]*v[:, i])
w, v = eig(a, left=1, right=0)
for i in range(3):
assert_array_almost_equal(a.T @ v[:, i], w[i]*v[:, i])
def test_simple_complex_eig(self):
a = array([[1, 2], [-2, 1]])
w, vl, vr = eig(a, left=1, right=1)
assert_array_almost_equal(w, array([1+2j, 1-2j]))
for i in range(2):
assert_array_almost_equal(a @ vr[:, i], w[i]*vr[:, i])
for i in range(2):
assert_array_almost_equal(a.conj().T @ vl[:, i],
w[i].conj()*vl[:, i])
def test_simple_complex(self):
a = array([[1, 2, 3], [1, 2, 3], [2, 5, 6+1j]])
w, vl, vr = eig(a, left=1, right=1)
for i in range(3):
assert_array_almost_equal(a @ vr[:, i], w[i]*vr[:, i])
for i in range(3):
assert_array_almost_equal(a.conj().T @ vl[:, i],
w[i].conj()*vl[:, i])
def test_gh_3054(self):
a = [[1]]
b = [[0]]
w, vr = eig(a, b, homogeneous_eigvals=True)
assert_allclose(w[1, 0], 0)
assert_(w[0, 0] != 0)
assert_allclose(vr, 1)
w, vr = eig(a, b)
assert_equal(w, np.inf)
assert_allclose(vr, 1)
def _check_gen_eig(self, A, B, atol_homog=1e-13, rtol_homog=1e-13):
if B is not None:
A, B = asarray(A), asarray(B)
B0 = B
else:
A = asarray(A)
B0 = B
B = np.eye(*A.shape)
msg = f"\n{A!r}\n{B!r}"
# Eigenvalues in homogeneous coordinates
w, vr = eig(A, B0, homogeneous_eigvals=True)
wt = eigvals(A, B0, homogeneous_eigvals=True)
val1 = A @ vr * w[1, :]
val2 = B @ vr * w[0, :]
for i in range(val1.shape[1]):
assert_allclose(val1[:, i], val2[:, i],
rtol=rtol_homog, atol=atol_homog, err_msg=msg)
if B0 is None:
assert_allclose(w[1, :], 1)
assert_allclose(wt[1, :], 1)
perm = np.lexsort(w)
permt = np.lexsort(wt)
assert_allclose(w[:, perm], wt[:, permt], atol=1e-7, rtol=1e-7,
err_msg=msg)
length = np.empty(len(vr))
for i in range(len(vr)):
length[i] = norm(vr[:, i])
assert_allclose(length, np.ones(length.size), err_msg=msg,
atol=1e-7, rtol=1e-7)
# Convert homogeneous coordinates
beta_nonzero = (w[1, :] != 0)
wh = w[0, beta_nonzero] / w[1, beta_nonzero]
# Eigenvalues in standard coordinates
w, vr = eig(A, B0)
wt = eigvals(A, B0)
val1 = A @ vr
val2 = B @ vr * w
res = val1 - val2
for i in range(res.shape[1]):
if np.all(isfinite(res[:, i])):
assert_allclose(res[:, i], 0,
rtol=1e-13, atol=1e-13, err_msg=msg)
# try to consistently order eigenvalues, including complex conjugate pairs
w_fin = w[isfinite(w)]
wt_fin = wt[isfinite(wt)]
# prune noise in the real parts
w_fin = -1j * np.real_if_close(1j*w_fin, tol=1e-10)
wt_fin = -1j * np.real_if_close(1j*wt_fin, tol=1e-10)
perm = argsort(abs(w_fin) + w_fin.imag)
permt = argsort(abs(wt_fin) + wt_fin.imag)
assert_allclose(w_fin[perm], wt_fin[permt],
atol=1e-7, rtol=1e-7, err_msg=msg)
length = np.empty(len(vr))
for i in range(len(vr)):
length[i] = norm(vr[:, i])
assert_allclose(length, np.ones(length.size), err_msg=msg)
# Compare homogeneous and nonhomogeneous versions
assert_allclose(sort(wh), sort(w[np.isfinite(w)]))
def test_singular(self):
# Example taken from
# https://web.archive.org/web/20040903121217/http://www.cs.umu.se/research/nla/singular_pairs/guptri/matlab.html
A = array([[22, 34, 31, 31, 17],
[45, 45, 42, 19, 29],
[39, 47, 49, 26, 34],
[27, 31, 26, 21, 15],
[38, 44, 44, 24, 30]])
B = array([[13, 26, 25, 17, 24],
[31, 46, 40, 26, 37],
[26, 40, 19, 25, 25],
[16, 25, 27, 14, 23],
[24, 35, 18, 21, 22]])
with np.errstate(all='ignore'):
self._check_gen_eig(A, B, atol_homog=5e-13)
def test_falker(self):
# Test matrices giving some Nan generalized eigenvalues.
M = diag(array([1, 0, 3]))
K = array(([2, -1, -1], [-1, 2, -1], [-1, -1, 2]))
D = array(([1, -1, 0], [-1, 1, 0], [0, 0, 0]))
Z = zeros((3, 3))
I3 = eye(3)
A = np.block([[I3, Z], [Z, -K]])
B = np.block([[Z, I3], [M, D]])
with np.errstate(all='ignore'):
self._check_gen_eig(A, B)
def test_bad_geneig(self):
# Ticket #709 (strange return values from DGGEV)
def matrices(omega):
c1 = -9 + omega**2
c2 = 2*omega
A = [[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, c1, 0],
[0, 0, 0, c1]]
B = [[0, 0, 1, 0],
[0, 0, 0, 1],
[1, 0, 0, -c2],
[0, 1, c2, 0]]
return A, B
# With a buggy LAPACK, this can fail for different omega on different
# machines -- so we need to test several values
with np.errstate(all='ignore'):
for k in range(100):
A, B = matrices(omega=k*5./100)
self._check_gen_eig(A, B)
def test_make_eigvals(self):
# Step through all paths in _make_eigvals
# Real eigenvalues
rng = np.random.RandomState(1234)
A = symrand(3, rng)
self._check_gen_eig(A, None)
B = symrand(3, rng)
self._check_gen_eig(A, B)
# Complex eigenvalues
A = rng.random((3, 3)) + 1j*rng.random((3, 3))
self._check_gen_eig(A, None)
B = rng.random((3, 3)) + 1j*rng.random((3, 3))
self._check_gen_eig(A, B)
def test_check_finite(self):
a = [[1, 2, 3], [1, 2, 3], [2, 5, 6]]
w, v = eig(a, check_finite=False)
exact_w = [(9+sqrt(93))/2, 0, (9-sqrt(93))/2]
v0 = array([1, 1, (1+sqrt(93)/3)/2])
v1 = array([3., 0, -1])
v2 = array([1, 1, (1-sqrt(93)/3)/2])
v0 = v0 / norm(v0)
v1 = v1 / norm(v1)
v2 = v2 / norm(v2)
assert_array_almost_equal(w, exact_w)
assert_array_almost_equal(v0, v[:, 0]*sign(v[0, 0]))
assert_array_almost_equal(v1, v[:, 1]*sign(v[0, 1]))
assert_array_almost_equal(v2, v[:, 2]*sign(v[0, 2]))
for i in range(3):
assert_array_almost_equal(a @ v[:, i], w[i]*v[:, i])
def test_not_square_error(self):
"""Check that passing a non-square array raises a ValueError."""
A = np.arange(6).reshape(3, 2)
assert_raises(ValueError, eig, A)
def test_shape_mismatch(self):
"""Check that passing arrays of with different shapes
raises a ValueError."""
A = eye(2)
B = np.arange(9.0).reshape(3, 3)
assert_raises(ValueError, eig, A, B)
assert_raises(ValueError, eig, B, A)
def test_gh_11577(self):
# https://github.com/scipy/scipy/issues/11577
# `A - lambda B` should have 4 and 8 among the eigenvalues, and this
# was apparently broken on some platforms
A = np.array([[12.0, 28.0, 76.0, 220.0],
[16.0, 32.0, 80.0, 224.0],
[24.0, 40.0, 88.0, 232.0],
[40.0, 56.0, 104.0, 248.0]], dtype='float64')
B = np.array([[2.0, 4.0, 10.0, 28.0],
[3.0, 5.0, 11.0, 29.0],
[5.0, 7.0, 13.0, 31.0],
[9.0, 11.0, 17.0, 35.0]], dtype='float64')
D, V = eig(A, B)
# The problem is ill-conditioned, and two other eigenvalues
# depend on ATLAS/OpenBLAS version, compiler version etc
# see gh-11577 for discussion
#
# NB: it is tempting to use `assert_allclose(D[:2], [4, 8])` instead but
# the ordering of eigenvalues also comes out different on different
# systems depending on who knows what.
with np.testing.suppress_warnings() as sup:
# isclose chokes on inf/nan values
sup.filter(RuntimeWarning, "invalid value encountered in multiply")
assert np.isclose(D, 4.0, atol=1e-14).any()
assert np.isclose(D, 8.0, atol=1e-14).any()
class TestEigBanded:
def setup_method(self):
self.create_bandmat()
def create_bandmat(self):
"""Create the full matrix `self.fullmat` and
the corresponding band matrix `self.bandmat`."""
N = 10
self.KL = 2 # number of subdiagonals (below the diagonal)
self.KU = 2 # number of superdiagonals (above the diagonal)
# symmetric band matrix
self.sym_mat = (diag(full(N, 1.0))
+ diag(full(N-1, -1.0), -1) + diag(full(N-1, -1.0), 1)
+ diag(full(N-2, -2.0), -2) + diag(full(N-2, -2.0), 2))
# hermitian band matrix
self.herm_mat = (diag(full(N, -1.0))
+ 1j*diag(full(N-1, 1.0), -1)
- 1j*diag(full(N-1, 1.0), 1)
+ diag(full(N-2, -2.0), -2)
+ diag(full(N-2, -2.0), 2))
# general real band matrix
self.real_mat = (diag(full(N, 1.0))
+ diag(full(N-1, -1.0), -1) + diag(full(N-1, -3.0), 1)
+ diag(full(N-2, 2.0), -2) + diag(full(N-2, -2.0), 2))
# general complex band matrix
self.comp_mat = (1j*diag(full(N, 1.0))
+ diag(full(N-1, -1.0), -1)
+ 1j*diag(full(N-1, -3.0), 1)
+ diag(full(N-2, 2.0), -2)
+ diag(full(N-2, -2.0), 2))
# Eigenvalues and -vectors from linalg.eig
ew, ev = linalg.eig(self.sym_mat)
ew = ew.real
args = argsort(ew)
self.w_sym_lin = ew[args]
self.evec_sym_lin = ev[:, args]
ew, ev = linalg.eig(self.herm_mat)
ew = ew.real
args = argsort(ew)
self.w_herm_lin = ew[args]
self.evec_herm_lin = ev[:, args]
# Extract upper bands from symmetric and hermitian band matrices
# (for use in dsbevd, dsbevx, zhbevd, zhbevx
# and their single precision versions)
LDAB = self.KU + 1
self.bandmat_sym = zeros((LDAB, N), dtype=float)
self.bandmat_herm = zeros((LDAB, N), dtype=complex)
for i in range(LDAB):
self.bandmat_sym[LDAB-i-1, i:N] = diag(self.sym_mat, i)
self.bandmat_herm[LDAB-i-1, i:N] = diag(self.herm_mat, i)
# Extract bands from general real and complex band matrix
# (for use in dgbtrf, dgbtrs and their single precision versions)
LDAB = 2*self.KL + self.KU + 1
self.bandmat_real = zeros((LDAB, N), dtype=float)
self.bandmat_real[2*self.KL, :] = diag(self.real_mat) # diagonal
for i in range(self.KL):
# superdiagonals
self.bandmat_real[2*self.KL-1-i, i+1:N] = diag(self.real_mat, i+1)
# subdiagonals
self.bandmat_real[2*self.KL+1+i, 0:N-1-i] = diag(self.real_mat,
-i-1)
self.bandmat_comp = zeros((LDAB, N), dtype=complex)
self.bandmat_comp[2*self.KL, :] = diag(self.comp_mat) # diagonal
for i in range(self.KL):
# superdiagonals
self.bandmat_comp[2*self.KL-1-i, i+1:N] = diag(self.comp_mat, i+1)
# subdiagonals
self.bandmat_comp[2*self.KL+1+i, 0:N-1-i] = diag(self.comp_mat,
-i-1)
# absolute value for linear equation system A*x = b
self.b = 1.0*arange(N)
self.bc = self.b * (1 + 1j)
#####################################################################
def test_dsbev(self):
"""Compare dsbev eigenvalues and eigenvectors with
the result of linalg.eig."""
w, evec, info = dsbev(self.bandmat_sym, compute_v=1)
evec_ = evec[:, argsort(w)]
assert_array_almost_equal(sort(w), self.w_sym_lin)
assert_array_almost_equal(abs(evec_), abs(self.evec_sym_lin))
def test_dsbevd(self):
"""Compare dsbevd eigenvalues and eigenvectors with
the result of linalg.eig."""
w, evec, info = dsbevd(self.bandmat_sym, compute_v=1)
evec_ = evec[:, argsort(w)]
assert_array_almost_equal(sort(w), self.w_sym_lin)
assert_array_almost_equal(abs(evec_), abs(self.evec_sym_lin))
def test_dsbevx(self):
"""Compare dsbevx eigenvalues and eigenvectors
with the result of linalg.eig."""
N, N = shape(self.sym_mat)
# Achtung: Argumente 0.0,0.0,range?
w, evec, num, ifail, info = dsbevx(self.bandmat_sym, 0.0, 0.0, 1, N,
compute_v=1, range=2)
evec_ = evec[:, argsort(w)]
assert_array_almost_equal(sort(w), self.w_sym_lin)
assert_array_almost_equal(abs(evec_), abs(self.evec_sym_lin))
def test_zhbevd(self):
"""Compare zhbevd eigenvalues and eigenvectors
with the result of linalg.eig."""
w, evec, info = zhbevd(self.bandmat_herm, compute_v=1)
evec_ = evec[:, argsort(w)]
assert_array_almost_equal(sort(w), self.w_herm_lin)
assert_array_almost_equal(abs(evec_), abs(self.evec_herm_lin))
def test_zhbevx(self):
"""Compare zhbevx eigenvalues and eigenvectors
with the result of linalg.eig."""
N, N = shape(self.herm_mat)
# Achtung: Argumente 0.0,0.0,range?
w, evec, num, ifail, info = zhbevx(self.bandmat_herm, 0.0, 0.0, 1, N,
compute_v=1, range=2)
evec_ = evec[:, argsort(w)]
assert_array_almost_equal(sort(w), self.w_herm_lin)
assert_array_almost_equal(abs(evec_), abs(self.evec_herm_lin))
def test_eigvals_banded(self):
"""Compare eigenvalues of eigvals_banded with those of linalg.eig."""
w_sym = eigvals_banded(self.bandmat_sym)
w_sym = w_sym.real
assert_array_almost_equal(sort(w_sym), self.w_sym_lin)
w_herm = eigvals_banded(self.bandmat_herm)
w_herm = w_herm.real
assert_array_almost_equal(sort(w_herm), self.w_herm_lin)
# extracting eigenvalues with respect to an index range
ind1 = 2
ind2 = np.longlong(6)
w_sym_ind = eigvals_banded(self.bandmat_sym,
select='i', select_range=(ind1, ind2))
assert_array_almost_equal(sort(w_sym_ind),
self.w_sym_lin[ind1:ind2+1])
w_herm_ind = eigvals_banded(self.bandmat_herm,
select='i', select_range=(ind1, ind2))
assert_array_almost_equal(sort(w_herm_ind),
self.w_herm_lin[ind1:ind2+1])
# extracting eigenvalues with respect to a value range
v_lower = self.w_sym_lin[ind1] - 1.0e-5
v_upper = self.w_sym_lin[ind2] + 1.0e-5
w_sym_val = eigvals_banded(self.bandmat_sym,
select='v', select_range=(v_lower, v_upper))
assert_array_almost_equal(sort(w_sym_val),
self.w_sym_lin[ind1:ind2+1])
v_lower = self.w_herm_lin[ind1] - 1.0e-5
v_upper = self.w_herm_lin[ind2] + 1.0e-5
w_herm_val = eigvals_banded(self.bandmat_herm,
select='v',
select_range=(v_lower, v_upper))
assert_array_almost_equal(sort(w_herm_val),
self.w_herm_lin[ind1:ind2+1])
w_sym = eigvals_banded(self.bandmat_sym, check_finite=False)
w_sym = w_sym.real
assert_array_almost_equal(sort(w_sym), self.w_sym_lin)
def test_eig_banded(self):
"""Compare eigenvalues and eigenvectors of eig_banded
with those of linalg.eig. """
w_sym, evec_sym = eig_banded(self.bandmat_sym)
evec_sym_ = evec_sym[:, argsort(w_sym.real)]
assert_array_almost_equal(sort(w_sym), self.w_sym_lin)
assert_array_almost_equal(abs(evec_sym_), abs(self.evec_sym_lin))
w_herm, evec_herm = eig_banded(self.bandmat_herm)
evec_herm_ = evec_herm[:, argsort(w_herm.real)]
assert_array_almost_equal(sort(w_herm), self.w_herm_lin)
assert_array_almost_equal(abs(evec_herm_), abs(self.evec_herm_lin))
# extracting eigenvalues with respect to an index range
ind1 = 2
ind2 = 6
w_sym_ind, evec_sym_ind = eig_banded(self.bandmat_sym,
select='i',
select_range=(ind1, ind2))
assert_array_almost_equal(sort(w_sym_ind),
self.w_sym_lin[ind1:ind2+1])
assert_array_almost_equal(abs(evec_sym_ind),
abs(self.evec_sym_lin[:, ind1:ind2+1]))
w_herm_ind, evec_herm_ind = eig_banded(self.bandmat_herm,
select='i',
select_range=(ind1, ind2))
assert_array_almost_equal(sort(w_herm_ind),
self.w_herm_lin[ind1:ind2+1])
assert_array_almost_equal(abs(evec_herm_ind),
abs(self.evec_herm_lin[:, ind1:ind2+1]))
# extracting eigenvalues with respect to a value range
v_lower = self.w_sym_lin[ind1] - 1.0e-5
v_upper = self.w_sym_lin[ind2] + 1.0e-5
w_sym_val, evec_sym_val = eig_banded(self.bandmat_sym,
select='v',
select_range=(v_lower, v_upper))
assert_array_almost_equal(sort(w_sym_val),
self.w_sym_lin[ind1:ind2+1])
assert_array_almost_equal(abs(evec_sym_val),
abs(self.evec_sym_lin[:, ind1:ind2+1]))
v_lower = self.w_herm_lin[ind1] - 1.0e-5
v_upper = self.w_herm_lin[ind2] + 1.0e-5
w_herm_val, evec_herm_val = eig_banded(self.bandmat_herm,
select='v',
select_range=(v_lower, v_upper))
assert_array_almost_equal(sort(w_herm_val),
self.w_herm_lin[ind1:ind2+1])
assert_array_almost_equal(abs(evec_herm_val),
abs(self.evec_herm_lin[:, ind1:ind2+1]))
w_sym, evec_sym = eig_banded(self.bandmat_sym, check_finite=False)
evec_sym_ = evec_sym[:, argsort(w_sym.real)]
assert_array_almost_equal(sort(w_sym), self.w_sym_lin)
assert_array_almost_equal(abs(evec_sym_), abs(self.evec_sym_lin))
def test_dgbtrf(self):
"""Compare dgbtrf LU factorisation with the LU factorisation result
of linalg.lu."""
M, N = shape(self.real_mat)
lu_symm_band, ipiv, info = dgbtrf(self.bandmat_real, self.KL, self.KU)
# extract matrix u from lu_symm_band
u = diag(lu_symm_band[2*self.KL, :])
for i in range(self.KL + self.KU):
u += diag(lu_symm_band[2*self.KL-1-i, i+1:N], i+1)
p_lin, l_lin, u_lin = lu(self.real_mat, permute_l=0)
assert_array_almost_equal(u, u_lin)
def test_zgbtrf(self):
"""Compare zgbtrf LU factorisation with the LU factorisation result
of linalg.lu."""
M, N = shape(self.comp_mat)
lu_symm_band, ipiv, info = zgbtrf(self.bandmat_comp, self.KL, self.KU)
# extract matrix u from lu_symm_band
u = diag(lu_symm_band[2*self.KL, :])
for i in range(self.KL + self.KU):
u += diag(lu_symm_band[2*self.KL-1-i, i+1:N], i+1)
p_lin, l_lin, u_lin = lu(self.comp_mat, permute_l=0)
assert_array_almost_equal(u, u_lin)
def test_dgbtrs(self):
"""Compare dgbtrs solutions for linear equation system A*x = b
with solutions of linalg.solve."""
lu_symm_band, ipiv, info = dgbtrf(self.bandmat_real, self.KL, self.KU)
y, info = dgbtrs(lu_symm_band, self.KL, self.KU, self.b, ipiv)
y_lin = linalg.solve(self.real_mat, self.b)
assert_array_almost_equal(y, y_lin)
def test_zgbtrs(self):
"""Compare zgbtrs solutions for linear equation system A*x = b
with solutions of linalg.solve."""
lu_symm_band, ipiv, info = zgbtrf(self.bandmat_comp, self.KL, self.KU)
y, info = zgbtrs(lu_symm_band, self.KL, self.KU, self.bc, ipiv)
y_lin = linalg.solve(self.comp_mat, self.bc)
assert_array_almost_equal(y, y_lin)
class TestEigTridiagonal:
def setup_method(self):
self.create_trimat()
def create_trimat(self):
"""Create the full matrix `self.fullmat`, `self.d`, and `self.e`."""
N = 10
# symmetric band matrix
self.d = full(N, 1.0)
self.e = full(N-1, -1.0)
self.full_mat = (diag(self.d) + diag(self.e, -1) + diag(self.e, 1))
ew, ev = linalg.eig(self.full_mat)
ew = ew.real
args = argsort(ew)
self.w = ew[args]
self.evec = ev[:, args]
def test_degenerate(self):
"""Test error conditions."""
# Wrong sizes
assert_raises(ValueError, eigvalsh_tridiagonal, self.d, self.e[:-1])
# Must be real
assert_raises(TypeError, eigvalsh_tridiagonal, self.d, self.e * 1j)
# Bad driver
assert_raises(TypeError, eigvalsh_tridiagonal, self.d, self.e,
lapack_driver=1.)
assert_raises(ValueError, eigvalsh_tridiagonal, self.d, self.e,
lapack_driver='foo')
# Bad bounds
assert_raises(ValueError, eigvalsh_tridiagonal, self.d, self.e,
select='i', select_range=(0, -1))
def test_eigvalsh_tridiagonal(self):
"""Compare eigenvalues of eigvalsh_tridiagonal with those of eig."""
# can't use ?STERF with subselection
for driver in ('sterf', 'stev', 'stebz', 'stemr', 'auto'):
w = eigvalsh_tridiagonal(self.d, self.e, lapack_driver=driver)
assert_array_almost_equal(sort(w), self.w)
for driver in ('sterf', 'stev'):
assert_raises(ValueError, eigvalsh_tridiagonal, self.d, self.e,
lapack_driver='stev', select='i',
select_range=(0, 1))
for driver in ('stebz', 'stemr', 'auto'):
# extracting eigenvalues with respect to the full index range
w_ind = eigvalsh_tridiagonal(
self.d, self.e, select='i', select_range=(0, len(self.d)-1),
lapack_driver=driver)
assert_array_almost_equal(sort(w_ind), self.w)
# extracting eigenvalues with respect to an index range
ind1 = 2
ind2 = 6
w_ind = eigvalsh_tridiagonal(
self.d, self.e, select='i', select_range=(ind1, ind2),
lapack_driver=driver)
assert_array_almost_equal(sort(w_ind), self.w[ind1:ind2+1])
# extracting eigenvalues with respect to a value range
v_lower = self.w[ind1] - 1.0e-5
v_upper = self.w[ind2] + 1.0e-5
w_val = eigvalsh_tridiagonal(
self.d, self.e, select='v', select_range=(v_lower, v_upper),
lapack_driver=driver)
assert_array_almost_equal(sort(w_val), self.w[ind1:ind2+1])
def test_eigh_tridiagonal(self):
"""Compare eigenvalues and eigenvectors of eigh_tridiagonal
with those of eig. """
# can't use ?STERF when eigenvectors are requested
assert_raises(ValueError, eigh_tridiagonal, self.d, self.e,
lapack_driver='sterf')
for driver in ('stebz', 'stev', 'stemr', 'auto'):
w, evec = eigh_tridiagonal(self.d, self.e, lapack_driver=driver)
evec_ = evec[:, argsort(w)]
assert_array_almost_equal(sort(w), self.w)
assert_array_almost_equal(abs(evec_), abs(self.evec))
assert_raises(ValueError, eigh_tridiagonal, self.d, self.e,
lapack_driver='stev', select='i', select_range=(0, 1))
for driver in ('stebz', 'stemr', 'auto'):
# extracting eigenvalues with respect to an index range
ind1 = 0
ind2 = len(self.d)-1
w, evec = eigh_tridiagonal(
self.d, self.e, select='i', select_range=(ind1, ind2),
lapack_driver=driver)
assert_array_almost_equal(sort(w), self.w)
assert_array_almost_equal(abs(evec), abs(self.evec))
ind1 = 2
ind2 = 6
w, evec = eigh_tridiagonal(
self.d, self.e, select='i', select_range=(ind1, ind2),
lapack_driver=driver)
assert_array_almost_equal(sort(w), self.w[ind1:ind2+1])
assert_array_almost_equal(abs(evec),
abs(self.evec[:, ind1:ind2+1]))
# extracting eigenvalues with respect to a value range
v_lower = self.w[ind1] - 1.0e-5
v_upper = self.w[ind2] + 1.0e-5
w, evec = eigh_tridiagonal(
self.d, self.e, select='v', select_range=(v_lower, v_upper),
lapack_driver=driver)
assert_array_almost_equal(sort(w), self.w[ind1:ind2+1])
assert_array_almost_equal(abs(evec),
abs(self.evec[:, ind1:ind2+1]))
def test_eigh_tridiagonal_1x1(self):
"""See gh-20075"""
a = np.array([-2.0])
b = np.array([])
x = eigh_tridiagonal(a, b, eigvals_only=True)
assert x.ndim == 1
assert_allclose(x, a)
x, V = eigh_tridiagonal(a, b, select="i", select_range=(0, 0))
assert x.ndim == 1
assert V.ndim == 2
assert_allclose(x, a)
assert_allclose(V, array([[1.]]))
x, V = eigh_tridiagonal(a, b, select="v", select_range=(-2, 0))
assert x.size == 0
assert x.shape == (0,)
assert V.shape == (1, 0)
class TestEigh:
def setup_class(self):
np.random.seed(1234)
def test_wrong_inputs(self):
# Nonsquare a
assert_raises(ValueError, eigh, np.ones([1, 2]))
# Nonsquare b
assert_raises(ValueError, eigh, np.ones([2, 2]), np.ones([2, 1]))
# Incompatible a, b sizes
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([2, 2]))
# Wrong type parameter for generalized problem
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
type=4)
# Both value and index subsets requested
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
subset_by_value=[1, 2], subset_by_index=[2, 4])
with np.testing.suppress_warnings() as sup:
sup.filter(DeprecationWarning, "Keyword argument 'eigvals")
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
subset_by_value=[1, 2], eigvals=[2, 4])
# Invalid upper index spec
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
subset_by_index=[0, 4])
with np.testing.suppress_warnings() as sup:
sup.filter(DeprecationWarning, "Keyword argument 'eigvals")
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
eigvals=[0, 4])
# Invalid lower index
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
subset_by_index=[-2, 2])
with np.testing.suppress_warnings() as sup:
sup.filter(DeprecationWarning, "Keyword argument 'eigvals")
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
eigvals=[-2, 2])
# Invalid index spec #2
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
subset_by_index=[2, 0])
with np.testing.suppress_warnings() as sup:
sup.filter(DeprecationWarning, "Keyword argument 'eigvals")
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
subset_by_index=[2, 0])
# Invalid value spec
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
subset_by_value=[2, 0])
# Invalid driver name
assert_raises(ValueError, eigh, np.ones([2, 2]), driver='wrong')
# Generalized driver selection without b
assert_raises(ValueError, eigh, np.ones([3, 3]), None, driver='gvx')
# Standard driver with b
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
driver='evr')
# Subset request from invalid driver
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
driver='gvd', subset_by_index=[1, 2])
assert_raises(ValueError, eigh, np.ones([3, 3]), np.ones([3, 3]),
driver='gvd', subset_by_index=[1, 2])
def test_nonpositive_b(self):
assert_raises(LinAlgError, eigh, np.ones([3, 3]), np.ones([3, 3]))
# index based subsets are done in the legacy test_eigh()
def test_value_subsets(self):
for ind, dt in enumerate(DTYPES):
a = _random_hermitian_matrix(20, dtype=dt)
w, v = eigh(a, subset_by_value=[-2, 2])
assert_equal(v.shape[1], len(w))
assert all((w > -2) & (w < 2))
b = _random_hermitian_matrix(20, posdef=True, dtype=dt)
w, v = eigh(a, b, subset_by_value=[-2, 2])
assert_equal(v.shape[1], len(w))
assert all((w > -2) & (w < 2))
def test_eigh_integer(self):
a = array([[1, 2], [2, 7]])
b = array([[3, 1], [1, 5]])
w, z = eigh(a)
w, z = eigh(a, b)
def test_eigh_of_sparse(self):
# This tests the rejection of inputs that eigh cannot currently handle.
import scipy.sparse
a = scipy.sparse.identity(2).tocsc()
b = np.atleast_2d(a)
assert_raises(ValueError, eigh, a)
assert_raises(ValueError, eigh, b)
@pytest.mark.parametrize('dtype_', DTYPES)
@pytest.mark.parametrize('driver', ("ev", "evd", "evr", "evx"))
def test_various_drivers_standard(self, driver, dtype_):
a = _random_hermitian_matrix(n=20, dtype=dtype_)
w, v = eigh(a, driver=driver)
assert_allclose(a @ v - (v * w), 0.,
atol=1000*np.finfo(dtype_).eps,
rtol=0.)
@pytest.mark.parametrize('driver', ("ev", "evd", "evr", "evx"))
def test_1x1_lwork(self, driver):
w, v = eigh([[1]], driver=driver)
assert_allclose(w, array([1.]), atol=1e-15)
assert_allclose(v, array([[1.]]), atol=1e-15)
# complex case now
w, v = eigh([[1j]], driver=driver)
assert_allclose(w, array([0]), atol=1e-15)
assert_allclose(v, array([[1.]]), atol=1e-15)
@pytest.mark.parametrize('type', (1, 2, 3))
@pytest.mark.parametrize('driver', ("gv", "gvd", "gvx"))
def test_various_drivers_generalized(self, driver, type):
atol = np.spacing(5000.)
a = _random_hermitian_matrix(20)
b = _random_hermitian_matrix(20, posdef=True)
w, v = eigh(a=a, b=b, driver=driver, type=type)
if type == 1:
assert_allclose(a @ v - w*(b @ v), 0., atol=atol, rtol=0.)
elif type == 2:
assert_allclose(a @ b @ v - v * w, 0., atol=atol, rtol=0.)
else:
assert_allclose(b @ a @ v - v * w, 0., atol=atol, rtol=0.)
def test_eigvalsh_new_args(self):
a = _random_hermitian_matrix(5)
w = eigvalsh(a, subset_by_index=[1, 2])
assert_equal(len(w), 2)
w2 = eigvalsh(a, subset_by_index=[1, 2])
assert_equal(len(w2), 2)
assert_allclose(w, w2)
b = np.diag([1, 1.2, 1.3, 1.5, 2])
w3 = eigvalsh(b, subset_by_value=[1, 1.4])
assert_equal(len(w3), 2)
assert_allclose(w3, np.array([1.2, 1.3]))
@pytest.mark.parametrize("method", [eigh, eigvalsh])
def test_deprecation_warnings(self, method):
with pytest.warns(DeprecationWarning,
match="Keyword argument 'turbo'"):
method(np.zeros((2, 2)), turbo=True)
with pytest.warns(DeprecationWarning,
match="Keyword argument 'eigvals'"):
method(np.zeros((2, 2)), eigvals=[0, 1])
with pytest.deprecated_call(match="use keyword arguments"):
method(np.zeros((2,2)), np.eye(2, 2), True)
def test_deprecation_results(self):
a = _random_hermitian_matrix(3)
b = _random_hermitian_matrix(3, posdef=True)
# check turbo gives same result as driver='gvd'
with np.testing.suppress_warnings() as sup:
sup.filter(DeprecationWarning, "Keyword argument 'turbo'")
w_dep, v_dep = eigh(a, b, turbo=True)
w, v = eigh(a, b, driver='gvd')
assert_allclose(w_dep, w)
assert_allclose(v_dep, v)
# check eigvals gives the same result as subset_by_index
with np.testing.suppress_warnings() as sup:
sup.filter(DeprecationWarning, "Keyword argument 'eigvals'")
w_dep, v_dep = eigh(a, eigvals=[0, 1])
w, v = eigh(a, subset_by_index=[0, 1])
assert_allclose(w_dep, w)
assert_allclose(v_dep, v)
class TestSVD_GESDD:
lapack_driver = 'gesdd'
def test_degenerate(self):
assert_raises(TypeError, svd, [[1.]], lapack_driver=1.)
assert_raises(ValueError, svd, [[1.]], lapack_driver='foo')
def test_simple(self):
a = [[1, 2, 3], [1, 20, 3], [2, 5, 6]]
for full_matrices in (True, False):
u, s, vh = svd(a, full_matrices=full_matrices,
lapack_driver=self.lapack_driver)
assert_array_almost_equal(u.T @ u, eye(3))
assert_array_almost_equal(vh.T @ vh, eye(3))
sigma = zeros((u.shape[0], vh.shape[0]), s.dtype.char)
for i in range(len(s)):
sigma[i, i] = s[i]
assert_array_almost_equal(u @ sigma @ vh, a)
def test_simple_singular(self):
a = [[1, 2, 3], [1, 2, 3], [2, 5, 6]]
for full_matrices in (True, False):
u, s, vh = svd(a, full_matrices=full_matrices,
lapack_driver=self.lapack_driver)
assert_array_almost_equal(u.T @ u, eye(3))
assert_array_almost_equal(vh.T @ vh, eye(3))
sigma = zeros((u.shape[0], vh.shape[0]), s.dtype.char)
for i in range(len(s)):
sigma[i, i] = s[i]
assert_array_almost_equal(u @ sigma @ vh, a)
def test_simple_underdet(self):
a = [[1, 2, 3], [4, 5, 6]]
for full_matrices in (True, False):
u, s, vh = svd(a, full_matrices=full_matrices,
lapack_driver=self.lapack_driver)
assert_array_almost_equal(u.T @ u, eye(u.shape[0]))
sigma = zeros((u.shape[0], vh.shape[0]), s.dtype.char)
for i in range(len(s)):
sigma[i, i] = s[i]
assert_array_almost_equal(u @ sigma @ vh, a)
def test_simple_overdet(self):
a = [[1, 2], [4, 5], [3, 4]]
for full_matrices in (True, False):
u, s, vh = svd(a, full_matrices=full_matrices,
lapack_driver=self.lapack_driver)
assert_array_almost_equal(u.T @ u, eye(u.shape[1]))
assert_array_almost_equal(vh.T @ vh, eye(2))
sigma = zeros((u.shape[1], vh.shape[0]), s.dtype.char)
for i in range(len(s)):
sigma[i, i] = s[i]
assert_array_almost_equal(u @ sigma @ vh, a)
def test_random(self):
rng = np.random.RandomState(1234)
n = 20
m = 15
for i in range(3):
for a in [rng.random([n, m]), rng.random([m, n])]:
for full_matrices in (True, False):
u, s, vh = svd(a, full_matrices=full_matrices,
lapack_driver=self.lapack_driver)
assert_array_almost_equal(u.T @ u, eye(u.shape[1]))
assert_array_almost_equal(vh @ vh.T, eye(vh.shape[0]))
sigma = zeros((u.shape[1], vh.shape[0]), s.dtype.char)
for i in range(len(s)):
sigma[i, i] = s[i]
assert_array_almost_equal(u @ sigma @ vh, a)
def test_simple_complex(self):
a = [[1, 2, 3], [1, 2j, 3], [2, 5, 6]]
for full_matrices in (True, False):
u, s, vh = svd(a, full_matrices=full_matrices,
lapack_driver=self.lapack_driver)
assert_array_almost_equal(u.conj().T @ u, eye(u.shape[1]))
assert_array_almost_equal(vh.conj().T @ vh, eye(vh.shape[0]))
sigma = zeros((u.shape[0], vh.shape[0]), s.dtype.char)
for i in range(len(s)):
sigma[i, i] = s[i]
assert_array_almost_equal(u @ sigma @ vh, a)
def test_random_complex(self):
rng = np.random.RandomState(1234)
n = 20
m = 15
for i in range(3):
for full_matrices in (True, False):
for a in [rng.random([n, m]), rng.random([m, n])]:
a = a + 1j*rng.random(list(a.shape))
u, s, vh = svd(a, full_matrices=full_matrices,
lapack_driver=self.lapack_driver)
assert_array_almost_equal(u.conj().T @ u,
eye(u.shape[1]))
# This fails when [m,n]
# assert_array_almost_equal(vh.conj().T @ vh,
# eye(len(vh),dtype=vh.dtype.char))
sigma = zeros((u.shape[1], vh.shape[0]), s.dtype.char)
for i in range(len(s)):
sigma[i, i] = s[i]
assert_array_almost_equal(u @ sigma @ vh, a)
def test_crash_1580(self):
rng = np.random.RandomState(1234)
sizes = [(13, 23), (30, 50), (60, 100)]
for sz in sizes:
for dt in [np.float32, np.float64, np.complex64, np.complex128]:
a = rng.rand(*sz).astype(dt)
# should not crash
svd(a, lapack_driver=self.lapack_driver)
def test_check_finite(self):
a = [[1, 2, 3], [1, 20, 3], [2, 5, 6]]
u, s, vh = svd(a, check_finite=False, lapack_driver=self.lapack_driver)
assert_array_almost_equal(u.T @ u, eye(3))
assert_array_almost_equal(vh.T @ vh, eye(3))
sigma = zeros((u.shape[0], vh.shape[0]), s.dtype.char)
for i in range(len(s)):
sigma[i, i] = s[i]
assert_array_almost_equal(u @ sigma @ vh, a)
def test_gh_5039(self):
# This is a smoke test for https://github.com/scipy/scipy/issues/5039
#
# The following is reported to raise "ValueError: On entry to DGESDD
# parameter number 12 had an illegal value".
# `interp1d([1,2,3,4], [1,2,3,4], kind='cubic')`
# This is reported to only show up on LAPACK 3.0.3.
#
# The matrix below is taken from the call to
# `B = _fitpack._bsplmat(order, xk)` in interpolate._find_smoothest
b = np.array(
[[0.16666667, 0.66666667, 0.16666667, 0., 0., 0.],
[0., 0.16666667, 0.66666667, 0.16666667, 0., 0.],
[0., 0., 0.16666667, 0.66666667, 0.16666667, 0.],
[0., 0., 0., 0.16666667, 0.66666667, 0.16666667]])
svd(b, lapack_driver=self.lapack_driver)
@pytest.mark.skipif(not HAS_ILP64, reason="64-bit LAPACK required")
@pytest.mark.slow
def test_large_matrix(self):
check_free_memory(free_mb=17000)
A = np.zeros([1, 2**31], dtype=np.float32)
A[0, -1] = 1
u, s, vh = svd(A, full_matrices=False)
assert_allclose(s[0], 1.0)
assert_allclose(u[0, 0] * vh[0, -1], 1.0)
class TestSVD_GESVD(TestSVD_GESDD):
lapack_driver = 'gesvd'
def test_svd_gesdd_nofegfault():
# svd(a) with {U,VT}.size > INT_MAX does not segfault
# cf https://github.com/scipy/scipy/issues/14001
df=np.ones((4799, 53130), dtype=np.float64)
with assert_raises(ValueError):
svd(df)
class TestSVDVals:
def test_empty(self):
for a in [[]], np.empty((2, 0)), np.ones((0, 3)):
s = svdvals(a)
assert_equal(s, np.empty(0))
def test_simple(self):
a = [[1, 2, 3], [1, 2, 3], [2, 5, 6]]
s = svdvals(a)
assert_(len(s) == 3)
assert_(s[0] >= s[1] >= s[2])
def test_simple_underdet(self):
a = [[1, 2, 3], [4, 5, 6]]
s = svdvals(a)
assert_(len(s) == 2)
assert_(s[0] >= s[1])
def test_simple_overdet(self):
a = [[1, 2], [4, 5], [3, 4]]
s = svdvals(a)
assert_(len(s) == 2)
assert_(s[0] >= s[1])
def test_simple_complex(self):
a = [[1, 2, 3], [1, 20, 3j], [2, 5, 6]]
s = svdvals(a)
assert_(len(s) == 3)
assert_(s[0] >= s[1] >= s[2])
def test_simple_underdet_complex(self):
a = [[1, 2, 3], [4, 5j, 6]]
s = svdvals(a)
assert_(len(s) == 2)
assert_(s[0] >= s[1])
def test_simple_overdet_complex(self):
a = [[1, 2], [4, 5], [3j, 4]]
s = svdvals(a)
assert_(len(s) == 2)
assert_(s[0] >= s[1])
def test_check_finite(self):
a = [[1, 2, 3], [1, 2, 3], [2, 5, 6]]
s = svdvals(a, check_finite=False)
assert_(len(s) == 3)
assert_(s[0] >= s[1] >= s[2])
@pytest.mark.slow
def test_crash_2609(self):
np.random.seed(1234)
a = np.random.rand(1500, 2800)
# Shouldn't crash:
svdvals(a)
class TestDiagSVD:
def test_simple(self):
assert_array_almost_equal(diagsvd([1, 0, 0], 3, 3),
[[1, 0, 0], [0, 0, 0], [0, 0, 0]])
class TestQR:
def test_simple(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
q, r = qr(a)
assert_array_almost_equal(q.T @ q, eye(3))
assert_array_almost_equal(q @ r, a)
def test_simple_left(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
q, r = qr(a)
c = [1, 2, 3]
qc, r2 = qr_multiply(a, c, "left")
assert_array_almost_equal(q @ c, qc)
assert_array_almost_equal(r, r2)
qc, r2 = qr_multiply(a, eye(3), "left")
assert_array_almost_equal(q, qc)
def test_simple_right(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
q, r = qr(a)
c = [1, 2, 3]
qc, r2 = qr_multiply(a, c)
assert_array_almost_equal(c @ q, qc)
assert_array_almost_equal(r, r2)
qc, r = qr_multiply(a, eye(3))
assert_array_almost_equal(q, qc)
def test_simple_pivoting(self):
a = np.asarray([[8, 2, 3], [2, 9, 3], [5, 3, 6]])
q, r, p = qr(a, pivoting=True)
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(3))
assert_array_almost_equal(q @ r, a[:, p])
q2, r2 = qr(a[:, p])
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_simple_left_pivoting(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
q, r, jpvt = qr(a, pivoting=True)
c = [1, 2, 3]
qc, r, jpvt = qr_multiply(a, c, "left", True)
assert_array_almost_equal(q @ c, qc)
def test_simple_right_pivoting(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
q, r, jpvt = qr(a, pivoting=True)
c = [1, 2, 3]
qc, r, jpvt = qr_multiply(a, c, pivoting=True)
assert_array_almost_equal(c @ q, qc)
def test_simple_trap(self):
a = [[8, 2, 3], [2, 9, 3]]
q, r = qr(a)
assert_array_almost_equal(q.T @ q, eye(2))
assert_array_almost_equal(q @ r, a)
def test_simple_trap_pivoting(self):
a = np.asarray([[8, 2, 3], [2, 9, 3]])
q, r, p = qr(a, pivoting=True)
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(2))
assert_array_almost_equal(q @ r, a[:, p])
q2, r2 = qr(a[:, p])
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_simple_tall(self):
# full version
a = [[8, 2], [2, 9], [5, 3]]
q, r = qr(a)
assert_array_almost_equal(q.T @ q, eye(3))
assert_array_almost_equal(q @ r, a)
def test_simple_tall_pivoting(self):
# full version pivoting
a = np.asarray([[8, 2], [2, 9], [5, 3]])
q, r, p = qr(a, pivoting=True)
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(3))
assert_array_almost_equal(q @ r, a[:, p])
q2, r2 = qr(a[:, p])
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_simple_tall_e(self):
# economy version
a = [[8, 2], [2, 9], [5, 3]]
q, r = qr(a, mode='economic')
assert_array_almost_equal(q.T @ q, eye(2))
assert_array_almost_equal(q @ r, a)
assert_equal(q.shape, (3, 2))
assert_equal(r.shape, (2, 2))
def test_simple_tall_e_pivoting(self):
# economy version pivoting
a = np.asarray([[8, 2], [2, 9], [5, 3]])
q, r, p = qr(a, pivoting=True, mode='economic')
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(2))
assert_array_almost_equal(q @ r, a[:, p])
q2, r2 = qr(a[:, p], mode='economic')
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_simple_tall_left(self):
a = [[8, 2], [2, 9], [5, 3]]
q, r = qr(a, mode="economic")
c = [1, 2]
qc, r2 = qr_multiply(a, c, "left")
assert_array_almost_equal(q @ c, qc)
assert_array_almost_equal(r, r2)
c = array([1, 2, 0])
qc, r2 = qr_multiply(a, c, "left", overwrite_c=True)
assert_array_almost_equal(q @ c[:2], qc)
qc, r = qr_multiply(a, eye(2), "left")
assert_array_almost_equal(qc, q)
def test_simple_tall_left_pivoting(self):
a = [[8, 2], [2, 9], [5, 3]]
q, r, jpvt = qr(a, mode="economic", pivoting=True)
c = [1, 2]
qc, r, kpvt = qr_multiply(a, c, "left", True)
assert_array_equal(jpvt, kpvt)
assert_array_almost_equal(q @ c, qc)
qc, r, jpvt = qr_multiply(a, eye(2), "left", True)
assert_array_almost_equal(qc, q)
def test_simple_tall_right(self):
a = [[8, 2], [2, 9], [5, 3]]
q, r = qr(a, mode="economic")
c = [1, 2, 3]
cq, r2 = qr_multiply(a, c)
assert_array_almost_equal(c @ q, cq)
assert_array_almost_equal(r, r2)
cq, r = qr_multiply(a, eye(3))
assert_array_almost_equal(cq, q)
def test_simple_tall_right_pivoting(self):
a = [[8, 2], [2, 9], [5, 3]]
q, r, jpvt = qr(a, pivoting=True, mode="economic")
c = [1, 2, 3]
cq, r, jpvt = qr_multiply(a, c, pivoting=True)
assert_array_almost_equal(c @ q, cq)
cq, r, jpvt = qr_multiply(a, eye(3), pivoting=True)
assert_array_almost_equal(cq, q)
def test_simple_fat(self):
# full version
a = [[8, 2, 5], [2, 9, 3]]
q, r = qr(a)
assert_array_almost_equal(q.T @ q, eye(2))
assert_array_almost_equal(q @ r, a)
assert_equal(q.shape, (2, 2))
assert_equal(r.shape, (2, 3))
def test_simple_fat_pivoting(self):
# full version pivoting
a = np.asarray([[8, 2, 5], [2, 9, 3]])
q, r, p = qr(a, pivoting=True)
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(2))
assert_array_almost_equal(q @ r, a[:, p])
assert_equal(q.shape, (2, 2))
assert_equal(r.shape, (2, 3))
q2, r2 = qr(a[:, p])
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_simple_fat_e(self):
# economy version
a = [[8, 2, 3], [2, 9, 5]]
q, r = qr(a, mode='economic')
assert_array_almost_equal(q.T @ q, eye(2))
assert_array_almost_equal(q @ r, a)
assert_equal(q.shape, (2, 2))
assert_equal(r.shape, (2, 3))
def test_simple_fat_e_pivoting(self):
# economy version pivoting
a = np.asarray([[8, 2, 3], [2, 9, 5]])
q, r, p = qr(a, pivoting=True, mode='economic')
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(2))
assert_array_almost_equal(q @ r, a[:, p])
assert_equal(q.shape, (2, 2))
assert_equal(r.shape, (2, 3))
q2, r2 = qr(a[:, p], mode='economic')
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_simple_fat_left(self):
a = [[8, 2, 3], [2, 9, 5]]
q, r = qr(a, mode="economic")
c = [1, 2]
qc, r2 = qr_multiply(a, c, "left")
assert_array_almost_equal(q @ c, qc)
assert_array_almost_equal(r, r2)
qc, r = qr_multiply(a, eye(2), "left")
assert_array_almost_equal(qc, q)
def test_simple_fat_left_pivoting(self):
a = [[8, 2, 3], [2, 9, 5]]
q, r, jpvt = qr(a, mode="economic", pivoting=True)
c = [1, 2]
qc, r, jpvt = qr_multiply(a, c, "left", True)
assert_array_almost_equal(q @ c, qc)
qc, r, jpvt = qr_multiply(a, eye(2), "left", True)
assert_array_almost_equal(qc, q)
def test_simple_fat_right(self):
a = [[8, 2, 3], [2, 9, 5]]
q, r = qr(a, mode="economic")
c = [1, 2]
cq, r2 = qr_multiply(a, c)
assert_array_almost_equal(c @ q, cq)
assert_array_almost_equal(r, r2)
cq, r = qr_multiply(a, eye(2))
assert_array_almost_equal(cq, q)
def test_simple_fat_right_pivoting(self):
a = [[8, 2, 3], [2, 9, 5]]
q, r, jpvt = qr(a, pivoting=True, mode="economic")
c = [1, 2]
cq, r, jpvt = qr_multiply(a, c, pivoting=True)
assert_array_almost_equal(c @ q, cq)
cq, r, jpvt = qr_multiply(a, eye(2), pivoting=True)
assert_array_almost_equal(cq, q)
def test_simple_complex(self):
a = [[3, 3+4j, 5], [5, 2, 2+7j], [3, 2, 7]]
q, r = qr(a)
assert_array_almost_equal(q.conj().T @ q, eye(3))
assert_array_almost_equal(q @ r, a)
def test_simple_complex_left(self):
a = [[3, 3+4j, 5], [5, 2, 2+7j], [3, 2, 7]]
q, r = qr(a)
c = [1, 2, 3+4j]
qc, r = qr_multiply(a, c, "left")
assert_array_almost_equal(q @ c, qc)
qc, r = qr_multiply(a, eye(3), "left")
assert_array_almost_equal(q, qc)
def test_simple_complex_right(self):
a = [[3, 3+4j, 5], [5, 2, 2+7j], [3, 2, 7]]
q, r = qr(a)
c = [1, 2, 3+4j]
qc, r = qr_multiply(a, c)
assert_array_almost_equal(c @ q, qc)
qc, r = qr_multiply(a, eye(3))
assert_array_almost_equal(q, qc)
def test_simple_tall_complex_left(self):
a = [[8, 2+3j], [2, 9], [5+7j, 3]]
q, r = qr(a, mode="economic")
c = [1, 2+2j]
qc, r2 = qr_multiply(a, c, "left")
assert_array_almost_equal(q @ c, qc)
assert_array_almost_equal(r, r2)
c = array([1, 2, 0])
qc, r2 = qr_multiply(a, c, "left", overwrite_c=True)
assert_array_almost_equal(q @ c[:2], qc)
qc, r = qr_multiply(a, eye(2), "left")
assert_array_almost_equal(qc, q)
def test_simple_complex_left_conjugate(self):
a = [[3, 3+4j, 5], [5, 2, 2+7j], [3, 2, 7]]
q, r = qr(a)
c = [1, 2, 3+4j]
qc, r = qr_multiply(a, c, "left", conjugate=True)
assert_array_almost_equal(q.conj() @ c, qc)
def test_simple_complex_tall_left_conjugate(self):
a = [[3, 3+4j], [5, 2+2j], [3, 2]]
q, r = qr(a, mode='economic')
c = [1, 3+4j]
qc, r = qr_multiply(a, c, "left", conjugate=True)
assert_array_almost_equal(q.conj() @ c, qc)
def test_simple_complex_right_conjugate(self):
a = [[3, 3+4j, 5], [5, 2, 2+7j], [3, 2, 7]]
q, r = qr(a)
c = np.array([1, 2, 3+4j])
qc, r = qr_multiply(a, c, conjugate=True)
assert_array_almost_equal(c @ q.conj(), qc)
def test_simple_complex_pivoting(self):
a = array([[3, 3+4j, 5], [5, 2, 2+7j], [3, 2, 7]])
q, r, p = qr(a, pivoting=True)
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.conj().T @ q, eye(3))
assert_array_almost_equal(q @ r, a[:, p])
q2, r2 = qr(a[:, p])
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_simple_complex_left_pivoting(self):
a = array([[3, 3+4j, 5], [5, 2, 2+7j], [3, 2, 7]])
q, r, jpvt = qr(a, pivoting=True)
c = [1, 2, 3+4j]
qc, r, jpvt = qr_multiply(a, c, "left", True)
assert_array_almost_equal(q @ c, qc)
def test_simple_complex_right_pivoting(self):
a = array([[3, 3+4j, 5], [5, 2, 2+7j], [3, 2, 7]])
q, r, jpvt = qr(a, pivoting=True)
c = [1, 2, 3+4j]
qc, r, jpvt = qr_multiply(a, c, pivoting=True)
assert_array_almost_equal(c @ q, qc)
def test_random(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n])
q, r = qr(a)
assert_array_almost_equal(q.T @ q, eye(n))
assert_array_almost_equal(q @ r, a)
def test_random_left(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n])
q, r = qr(a)
c = rng.random([n])
qc, r = qr_multiply(a, c, "left")
assert_array_almost_equal(q @ c, qc)
qc, r = qr_multiply(a, eye(n), "left")
assert_array_almost_equal(q, qc)
def test_random_right(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n])
q, r = qr(a)
c = rng.random([n])
cq, r = qr_multiply(a, c)
assert_array_almost_equal(c @ q, cq)
cq, r = qr_multiply(a, eye(n))
assert_array_almost_equal(q, cq)
def test_random_pivoting(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n])
q, r, p = qr(a, pivoting=True)
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(n))
assert_array_almost_equal(q @ r, a[:, p])
q2, r2 = qr(a[:, p])
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_random_tall(self):
rng = np.random.RandomState(1234)
# full version
m = 200
n = 100
for k in range(2):
a = rng.random([m, n])
q, r = qr(a)
assert_array_almost_equal(q.T @ q, eye(m))
assert_array_almost_equal(q @ r, a)
def test_random_tall_left(self):
rng = np.random.RandomState(1234)
# full version
m = 200
n = 100
for k in range(2):
a = rng.random([m, n])
q, r = qr(a, mode="economic")
c = rng.random([n])
qc, r = qr_multiply(a, c, "left")
assert_array_almost_equal(q @ c, qc)
qc, r = qr_multiply(a, eye(n), "left")
assert_array_almost_equal(qc, q)
def test_random_tall_right(self):
rng = np.random.RandomState(1234)
# full version
m = 200
n = 100
for k in range(2):
a = rng.random([m, n])
q, r = qr(a, mode="economic")
c = rng.random([m])
cq, r = qr_multiply(a, c)
assert_array_almost_equal(c @ q, cq)
cq, r = qr_multiply(a, eye(m))
assert_array_almost_equal(cq, q)
def test_random_tall_pivoting(self):
rng = np.random.RandomState(1234)
# full version pivoting
m = 200
n = 100
for k in range(2):
a = rng.random([m, n])
q, r, p = qr(a, pivoting=True)
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(m))
assert_array_almost_equal(q @ r, a[:, p])
q2, r2 = qr(a[:, p])
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_random_tall_e(self):
rng = np.random.RandomState(1234)
# economy version
m = 200
n = 100
for k in range(2):
a = rng.random([m, n])
q, r = qr(a, mode='economic')
assert_array_almost_equal(q.T @ q, eye(n))
assert_array_almost_equal(q @ r, a)
assert_equal(q.shape, (m, n))
assert_equal(r.shape, (n, n))
def test_random_tall_e_pivoting(self):
rng = np.random.RandomState(1234)
# economy version pivoting
m = 200
n = 100
for k in range(2):
a = rng.random([m, n])
q, r, p = qr(a, pivoting=True, mode='economic')
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(n))
assert_array_almost_equal(q @ r, a[:, p])
assert_equal(q.shape, (m, n))
assert_equal(r.shape, (n, n))
q2, r2 = qr(a[:, p], mode='economic')
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_random_trap(self):
rng = np.random.RandomState(1234)
m = 100
n = 200
for k in range(2):
a = rng.random([m, n])
q, r = qr(a)
assert_array_almost_equal(q.T @ q, eye(m))
assert_array_almost_equal(q @ r, a)
def test_random_trap_pivoting(self):
rng = np.random.RandomState(1234)
m = 100
n = 200
for k in range(2):
a = rng.random([m, n])
q, r, p = qr(a, pivoting=True)
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.T @ q, eye(m))
assert_array_almost_equal(q @ r, a[:, p])
q2, r2 = qr(a[:, p])
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_random_complex(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n]) + 1j*rng.random([n, n])
q, r = qr(a)
assert_array_almost_equal(q.conj().T @ q, eye(n))
assert_array_almost_equal(q @ r, a)
def test_random_complex_left(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n]) + 1j*rng.random([n, n])
q, r = qr(a)
c = rng.random([n]) + 1j*rng.random([n])
qc, r = qr_multiply(a, c, "left")
assert_array_almost_equal(q @ c, qc)
qc, r = qr_multiply(a, eye(n), "left")
assert_array_almost_equal(q, qc)
def test_random_complex_right(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n]) + 1j*rng.random([n, n])
q, r = qr(a)
c = rng.random([n]) + 1j*rng.random([n])
cq, r = qr_multiply(a, c)
assert_array_almost_equal(c @ q, cq)
cq, r = qr_multiply(a, eye(n))
assert_array_almost_equal(q, cq)
def test_random_complex_pivoting(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n]) + 1j*rng.random([n, n])
q, r, p = qr(a, pivoting=True)
d = abs(diag(r))
assert_(np.all(d[1:] <= d[:-1]))
assert_array_almost_equal(q.conj().T @ q, eye(n))
assert_array_almost_equal(q @ r, a[:, p])
q2, r2 = qr(a[:, p])
assert_array_almost_equal(q, q2)
assert_array_almost_equal(r, r2)
def test_check_finite(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
q, r = qr(a, check_finite=False)
assert_array_almost_equal(q.T @ q, eye(3))
assert_array_almost_equal(q @ r, a)
def test_lwork(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
# Get comparison values
q, r = qr(a, lwork=None)
# Test against minimum valid lwork
q2, r2 = qr(a, lwork=3)
assert_array_almost_equal(q2, q)
assert_array_almost_equal(r2, r)
# Test against larger lwork
q3, r3 = qr(a, lwork=10)
assert_array_almost_equal(q3, q)
assert_array_almost_equal(r3, r)
# Test against explicit lwork=-1
q4, r4 = qr(a, lwork=-1)
assert_array_almost_equal(q4, q)
assert_array_almost_equal(r4, r)
# Test against invalid lwork
assert_raises(Exception, qr, (a,), {'lwork': 0})
assert_raises(Exception, qr, (a,), {'lwork': 2})
class TestRQ:
def test_simple(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
r, q = rq(a)
assert_array_almost_equal(q @ q.T, eye(3))
assert_array_almost_equal(r @ q, a)
def test_r(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
r, q = rq(a)
r2 = rq(a, mode='r')
assert_array_almost_equal(r, r2)
def test_random(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n])
r, q = rq(a)
assert_array_almost_equal(q @ q.T, eye(n))
assert_array_almost_equal(r @ q, a)
def test_simple_trap(self):
a = [[8, 2, 3], [2, 9, 3]]
r, q = rq(a)
assert_array_almost_equal(q.T @ q, eye(3))
assert_array_almost_equal(r @ q, a)
def test_simple_tall(self):
a = [[8, 2], [2, 9], [5, 3]]
r, q = rq(a)
assert_array_almost_equal(q.T @ q, eye(2))
assert_array_almost_equal(r @ q, a)
def test_simple_fat(self):
a = [[8, 2, 5], [2, 9, 3]]
r, q = rq(a)
assert_array_almost_equal(q @ q.T, eye(3))
assert_array_almost_equal(r @ q, a)
def test_simple_complex(self):
a = [[3, 3+4j, 5], [5, 2, 2+7j], [3, 2, 7]]
r, q = rq(a)
assert_array_almost_equal(q @ q.conj().T, eye(3))
assert_array_almost_equal(r @ q, a)
def test_random_tall(self):
rng = np.random.RandomState(1234)
m = 200
n = 100
for k in range(2):
a = rng.random([m, n])
r, q = rq(a)
assert_array_almost_equal(q @ q.T, eye(n))
assert_array_almost_equal(r @ q, a)
def test_random_trap(self):
rng = np.random.RandomState(1234)
m = 100
n = 200
for k in range(2):
a = rng.random([m, n])
r, q = rq(a)
assert_array_almost_equal(q @ q.T, eye(n))
assert_array_almost_equal(r @ q, a)
def test_random_trap_economic(self):
rng = np.random.RandomState(1234)
m = 100
n = 200
for k in range(2):
a = rng.random([m, n])
r, q = rq(a, mode='economic')
assert_array_almost_equal(q @ q.T, eye(m))
assert_array_almost_equal(r @ q, a)
assert_equal(q.shape, (m, n))
assert_equal(r.shape, (m, m))
def test_random_complex(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n]) + 1j*rng.random([n, n])
r, q = rq(a)
assert_array_almost_equal(q @ q.conj().T, eye(n))
assert_array_almost_equal(r @ q, a)
def test_random_complex_economic(self):
rng = np.random.RandomState(1234)
m = 100
n = 200
for k in range(2):
a = rng.random([m, n]) + 1j*rng.random([m, n])
r, q = rq(a, mode='economic')
assert_array_almost_equal(q @ q.conj().T, eye(m))
assert_array_almost_equal(r @ q, a)
assert_equal(q.shape, (m, n))
assert_equal(r.shape, (m, m))
def test_check_finite(self):
a = [[8, 2, 3], [2, 9, 3], [5, 3, 6]]
r, q = rq(a, check_finite=False)
assert_array_almost_equal(q @ q.T, eye(3))
assert_array_almost_equal(r @ q, a)
class TestSchur:
def check_schur(self, a, t, u, rtol, atol):
# Check that the Schur decomposition is correct.
assert_allclose(u @ t @ u.conj().T, a, rtol=rtol, atol=atol,
err_msg="Schur decomposition does not match 'a'")
# The expected value of u @ u.H - I is all zeros, so test
# with absolute tolerance only.
assert_allclose(u @ u.conj().T - np.eye(len(u)), 0, rtol=0, atol=atol,
err_msg="u is not unitary")
def test_simple(self):
a = [[8, 12, 3], [2, 9, 3], [10, 3, 6]]
t, z = schur(a)
self.check_schur(a, t, z, rtol=1e-14, atol=5e-15)
tc, zc = schur(a, 'complex')
assert_(np.any(ravel(iscomplex(zc))) and np.any(ravel(iscomplex(tc))))
self.check_schur(a, tc, zc, rtol=1e-14, atol=5e-15)
tc2, zc2 = rsf2csf(tc, zc)
self.check_schur(a, tc2, zc2, rtol=1e-14, atol=5e-15)
@pytest.mark.parametrize(
'sort, expected_diag',
[('lhp', [-np.sqrt(2), -0.5, np.sqrt(2), 0.5]),
('rhp', [np.sqrt(2), 0.5, -np.sqrt(2), -0.5]),
('iuc', [-0.5, 0.5, np.sqrt(2), -np.sqrt(2)]),
('ouc', [np.sqrt(2), -np.sqrt(2), -0.5, 0.5]),
(lambda x: x >= 0.0, [np.sqrt(2), 0.5, -np.sqrt(2), -0.5])]
)
def test_sort(self, sort, expected_diag):
# The exact eigenvalues of this matrix are
# -sqrt(2), sqrt(2), -1/2, 1/2.
a = [[4., 3., 1., -1.],
[-4.5, -3.5, -1., 1.],
[9., 6., -4., 4.5],
[6., 4., -3., 3.5]]
t, u, sdim = schur(a, sort=sort)
self.check_schur(a, t, u, rtol=1e-14, atol=5e-15)
assert_allclose(np.diag(t), expected_diag, rtol=1e-12)
assert_equal(2, sdim)
def test_sort_errors(self):
a = [[4., 3., 1., -1.],
[-4.5, -3.5, -1., 1.],
[9., 6., -4., 4.5],
[6., 4., -3., 3.5]]
assert_raises(ValueError, schur, a, sort='unsupported')
assert_raises(ValueError, schur, a, sort=1)
def test_check_finite(self):
a = [[8, 12, 3], [2, 9, 3], [10, 3, 6]]
t, z = schur(a, check_finite=False)
assert_array_almost_equal(z @ t @ z.conj().T, a)
class TestHessenberg:
def test_simple(self):
a = [[-149, -50, -154],
[537, 180, 546],
[-27, -9, -25]]
h1 = [[-149.0000, 42.2037, -156.3165],
[-537.6783, 152.5511, -554.9272],
[0, 0.0728, 2.4489]]
h, q = hessenberg(a, calc_q=1)
assert_array_almost_equal(q.T @ a @ q, h)
assert_array_almost_equal(h, h1, decimal=4)
def test_simple_complex(self):
a = [[-149, -50, -154],
[537, 180j, 546],
[-27j, -9, -25]]
h, q = hessenberg(a, calc_q=1)
assert_array_almost_equal(q.conj().T @ a @ q, h)
def test_simple2(self):
a = [[1, 2, 3, 4, 5, 6, 7],
[0, 2, 3, 4, 6, 7, 2],
[0, 2, 2, 3, 0, 3, 2],
[0, 0, 2, 8, 0, 0, 2],
[0, 3, 1, 2, 0, 1, 2],
[0, 1, 2, 3, 0, 1, 0],
[0, 0, 0, 0, 0, 1, 2]]
h, q = hessenberg(a, calc_q=1)
assert_array_almost_equal(q.T @ a @ q, h)
def test_simple3(self):
a = np.eye(3)
a[-1, 0] = 2
h, q = hessenberg(a, calc_q=1)
assert_array_almost_equal(q.T @ a @ q, h)
def test_random(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n])
h, q = hessenberg(a, calc_q=1)
assert_array_almost_equal(q.T @ a @ q, h)
def test_random_complex(self):
rng = np.random.RandomState(1234)
n = 20
for k in range(2):
a = rng.random([n, n]) + 1j*rng.random([n, n])
h, q = hessenberg(a, calc_q=1)
assert_array_almost_equal(q.conj().T @ a @ q, h)
def test_check_finite(self):
a = [[-149, -50, -154],
[537, 180, 546],
[-27, -9, -25]]
h1 = [[-149.0000, 42.2037, -156.3165],
[-537.6783, 152.5511, -554.9272],
[0, 0.0728, 2.4489]]
h, q = hessenberg(a, calc_q=1, check_finite=False)
assert_array_almost_equal(q.T @ a @ q, h)
assert_array_almost_equal(h, h1, decimal=4)
def test_2x2(self):
a = [[2, 1], [7, 12]]
h, q = hessenberg(a, calc_q=1)
assert_array_almost_equal(q, np.eye(2))
assert_array_almost_equal(h, a)
b = [[2-7j, 1+2j], [7+3j, 12-2j]]
h2, q2 = hessenberg(b, calc_q=1)
assert_array_almost_equal(q2, np.eye(2))
assert_array_almost_equal(h2, b)
blas_provider = blas_version = None
if CONFIG is not None:
blas_provider = CONFIG['Build Dependencies']['blas']['name']
blas_version = CONFIG['Build Dependencies']['blas']['version']
class TestQZ:
@pytest.mark.xfail(
sys.platform == 'darwin' and
blas_provider == 'openblas' and
blas_version < "0.3.21.dev",
reason="gges[float32] broken for OpenBLAS on macOS, see gh-16949"
)
def test_qz_single(self):
rng = np.random.RandomState(12345)
n = 5
A = rng.random([n, n]).astype(float32)
B = rng.random([n, n]).astype(float32)
AA, BB, Q, Z = qz(A, B)
assert_array_almost_equal(Q @ AA @ Z.T, A, decimal=5)
assert_array_almost_equal(Q @ BB @ Z.T, B, decimal=5)
assert_array_almost_equal(Q @ Q.T, eye(n), decimal=5)
assert_array_almost_equal(Z @ Z.T, eye(n), decimal=5)
assert_(np.all(diag(BB) >= 0))
def test_qz_double(self):
rng = np.random.RandomState(12345)
n = 5
A = rng.random([n, n])
B = rng.random([n, n])
AA, BB, Q, Z = qz(A, B)
assert_array_almost_equal(Q @ AA @ Z.T, A)
assert_array_almost_equal(Q @ BB @ Z.T, B)
assert_array_almost_equal(Q @ Q.T, eye(n))
assert_array_almost_equal(Z @ Z.T, eye(n))
assert_(np.all(diag(BB) >= 0))
def test_qz_complex(self):
rng = np.random.RandomState(12345)
n = 5
A = rng.random([n, n]) + 1j*rng.random([n, n])
B = rng.random([n, n]) + 1j*rng.random([n, n])
AA, BB, Q, Z = qz(A, B)
assert_array_almost_equal(Q @ AA @ Z.conj().T, A)
assert_array_almost_equal(Q @ BB @ Z.conj().T, B)
assert_array_almost_equal(Q @ Q.conj().T, eye(n))
assert_array_almost_equal(Z @ Z.conj().T, eye(n))
assert_(np.all(diag(BB) >= 0))
assert_(np.all(diag(BB).imag == 0))
def test_qz_complex64(self):
rng = np.random.RandomState(12345)
n = 5
A = (rng.random([n, n]) + 1j*rng.random([n, n])).astype(complex64)
B = (rng.random([n, n]) + 1j*rng.random([n, n])).astype(complex64)
AA, BB, Q, Z = qz(A, B)
assert_array_almost_equal(Q @ AA @ Z.conj().T, A, decimal=5)
assert_array_almost_equal(Q @ BB @ Z.conj().T, B, decimal=5)
assert_array_almost_equal(Q @ Q.conj().T, eye(n), decimal=5)
assert_array_almost_equal(Z @ Z.conj().T, eye(n), decimal=5)
assert_(np.all(diag(BB) >= 0))
assert_(np.all(diag(BB).imag == 0))
def test_qz_double_complex(self):
rng = np.random.RandomState(12345)
n = 5
A = rng.random([n, n])
B = rng.random([n, n])
AA, BB, Q, Z = qz(A, B, output='complex')
aa = Q @ AA @ Z.conj().T
assert_array_almost_equal(aa.real, A)
assert_array_almost_equal(aa.imag, 0)
bb = Q @ BB @ Z.conj().T
assert_array_almost_equal(bb.real, B)
assert_array_almost_equal(bb.imag, 0)
assert_array_almost_equal(Q @ Q.conj().T, eye(n))
assert_array_almost_equal(Z @ Z.conj().T, eye(n))
assert_(np.all(diag(BB) >= 0))
def test_qz_double_sort(self):
# from https://www.nag.com/lapack-ex/node119.html
# NOTE: These matrices may be ill-conditioned and lead to a
# seg fault on certain python versions when compiled with
# sse2 or sse3 older ATLAS/LAPACK binaries for windows
# A = np.array([[3.9, 12.5, -34.5, -0.5],
# [ 4.3, 21.5, -47.5, 7.5],
# [ 4.3, 21.5, -43.5, 3.5],
# [ 4.4, 26.0, -46.0, 6.0 ]])
# B = np.array([[ 1.0, 2.0, -3.0, 1.0],
# [1.0, 3.0, -5.0, 4.0],
# [1.0, 3.0, -4.0, 3.0],
# [1.0, 3.0, -4.0, 4.0]])
A = np.array([[3.9, 12.5, -34.5, 2.5],
[4.3, 21.5, -47.5, 7.5],
[4.3, 1.5, -43.5, 3.5],
[4.4, 6.0, -46.0, 6.0]])
B = np.array([[1.0, 1.0, -3.0, 1.0],
[1.0, 3.0, -5.0, 4.4],
[1.0, 2.0, -4.0, 1.0],
[1.2, 3.0, -4.0, 4.0]])
assert_raises(ValueError, qz, A, B, sort=lambda ar, ai, beta: ai == 0)
if False:
AA, BB, Q, Z, sdim = qz(A, B, sort=lambda ar, ai, beta: ai == 0)
# assert_(sdim == 2)
assert_(sdim == 4)
assert_array_almost_equal(Q @ AA @ Z.T, A)
assert_array_almost_equal(Q @ BB @ Z.T, B)
# test absolute values bc the sign is ambiguous and
# might be platform dependent
assert_array_almost_equal(np.abs(AA), np.abs(np.array(
[[35.7864, -80.9061, -12.0629, -9.498],
[0., 2.7638, -2.3505, 7.3256],
[0., 0., 0.6258, -0.0398],
[0., 0., 0., -12.8217]])), 4)
assert_array_almost_equal(np.abs(BB), np.abs(np.array(
[[4.5324, -8.7878, 3.2357, -3.5526],
[0., 1.4314, -2.1894, 0.9709],
[0., 0., 1.3126, -0.3468],
[0., 0., 0., 0.559]])), 4)
assert_array_almost_equal(np.abs(Q), np.abs(np.array(
[[-0.4193, -0.605, -0.1894, -0.6498],
[-0.5495, 0.6987, 0.2654, -0.3734],
[-0.4973, -0.3682, 0.6194, 0.4832],
[-0.5243, 0.1008, -0.7142, 0.4526]])), 4)
assert_array_almost_equal(np.abs(Z), np.abs(np.array(
[[-0.9471, -0.2971, -0.1217, 0.0055],
[-0.0367, 0.1209, 0.0358, 0.9913],
[0.3171, -0.9041, -0.2547, 0.1312],
[0.0346, 0.2824, -0.9587, 0.0014]])), 4)
# test absolute values bc the sign is ambiguous and might be platform
# dependent
# assert_array_almost_equal(abs(AA), abs(np.array([
# [3.8009, -69.4505, 50.3135, -43.2884],
# [0.0000, 9.2033, -0.2001, 5.9881],
# [0.0000, 0.0000, 1.4279, 4.4453],
# [0.0000, 0.0000, 0.9019, -1.1962]])), 4)
# assert_array_almost_equal(abs(BB), abs(np.array([
# [1.9005, -10.2285, 0.8658, -5.2134],
# [0.0000, 2.3008, 0.7915, 0.4262],
# [0.0000, 0.0000, 0.8101, 0.0000],
# [0.0000, 0.0000, 0.0000, -0.2823]])), 4)
# assert_array_almost_equal(abs(Q), abs(np.array([
# [0.4642, 0.7886, 0.2915, -0.2786],
# [0.5002, -0.5986, 0.5638, -0.2713],
# [0.5002, 0.0154, -0.0107, 0.8657],
# [0.5331, -0.1395, -0.7727, -0.3151]])), 4)
# assert_array_almost_equal(dot(Q,Q.T), eye(4))
# assert_array_almost_equal(abs(Z), abs(np.array([
# [0.9961, -0.0014, 0.0887, -0.0026],
# [0.0057, -0.0404, -0.0938, -0.9948],
# [0.0626, 0.7194, -0.6908, 0.0363],
# [0.0626, -0.6934, -0.7114, 0.0956]])), 4)
# assert_array_almost_equal(dot(Z,Z.T), eye(4))
# def test_qz_complex_sort(self):
# cA = np.array([
# [-21.10+22.50*1j, 53.50+-50.50*1j, -34.50+127.50*1j, 7.50+ 0.50*1j],
# [-0.46+ -7.78*1j, -3.50+-37.50*1j, -15.50+ 58.50*1j,-10.50+ -1.50*1j],
# [ 4.30+ -5.50*1j, 39.70+-17.10*1j, -68.50+ 12.50*1j, -7.50+ -3.50*1j],
# [ 5.50+ 4.40*1j, 14.40+ 43.30*1j, -32.50+-46.00*1j,-19.00+-32.50*1j]])
# cB = np.array([
# [1.00+ -5.00*1j, 1.60+ 1.20*1j,-3.00+ 0.00*1j, 0.00+ -1.00*1j],
# [0.80+ -0.60*1j, 3.00+ -5.00*1j,-4.00+ 3.00*1j,-2.40+ -3.20*1j],
# [1.00+ 0.00*1j, 2.40+ 1.80*1j,-4.00+ -5.00*1j, 0.00+ -3.00*1j],
# [0.00+ 1.00*1j,-1.80+ 2.40*1j, 0.00+ -4.00*1j, 4.00+ -5.00*1j]])
# AAS,BBS,QS,ZS,sdim = qz(cA,cB,sort='lhp')
# eigenvalues = diag(AAS)/diag(BBS)
# assert_(np.all(np.real(eigenvalues[:sdim] < 0)))
# assert_(np.all(np.real(eigenvalues[sdim:] > 0)))
def test_check_finite(self):
rng = np.random.RandomState(12345)
n = 5
A = rng.random([n, n])
B = rng.random([n, n])
AA, BB, Q, Z = qz(A, B, check_finite=False)
assert_array_almost_equal(Q @ AA @ Z.T, A)
assert_array_almost_equal(Q @ BB @ Z.T, B)
assert_array_almost_equal(Q @ Q.T, eye(n))
assert_array_almost_equal(Z @ Z.T, eye(n))
assert_(np.all(diag(BB) >= 0))
class TestOrdQZ:
@classmethod
def setup_class(cls):
# https://www.nag.com/lapack-ex/node119.html
A1 = np.array([[-21.10 - 22.50j, 53.5 - 50.5j, -34.5 + 127.5j,
7.5 + 0.5j],
[-0.46 - 7.78j, -3.5 - 37.5j, -15.5 + 58.5j,
-10.5 - 1.5j],
[4.30 - 5.50j, 39.7 - 17.1j, -68.5 + 12.5j,
-7.5 - 3.5j],
[5.50 + 4.40j, 14.4 + 43.3j, -32.5 - 46.0j,
-19.0 - 32.5j]])
B1 = np.array([[1.0 - 5.0j, 1.6 + 1.2j, -3 + 0j, 0.0 - 1.0j],
[0.8 - 0.6j, .0 - 5.0j, -4 + 3j, -2.4 - 3.2j],
[1.0 + 0.0j, 2.4 + 1.8j, -4 - 5j, 0.0 - 3.0j],
[0.0 + 1.0j, -1.8 + 2.4j, 0 - 4j, 4.0 - 5.0j]])
# https://www.nag.com/numeric/fl/nagdoc_fl23/xhtml/F08/f08yuf.xml
A2 = np.array([[3.9, 12.5, -34.5, -0.5],
[4.3, 21.5, -47.5, 7.5],
[4.3, 21.5, -43.5, 3.5],
[4.4, 26.0, -46.0, 6.0]])
B2 = np.array([[1, 2, -3, 1],
[1, 3, -5, 4],
[1, 3, -4, 3],
[1, 3, -4, 4]])
# example with the eigenvalues
# -0.33891648, 1.61217396+0.74013521j, 1.61217396-0.74013521j,
# 0.61244091
# thus featuring:
# * one complex conjugate eigenvalue pair,
# * one eigenvalue in the lhp
# * 2 eigenvalues in the unit circle
# * 2 non-real eigenvalues
A3 = np.array([[5., 1., 3., 3.],
[4., 4., 2., 7.],
[7., 4., 1., 3.],
[0., 4., 8., 7.]])
B3 = np.array([[8., 10., 6., 10.],
[7., 7., 2., 9.],
[9., 1., 6., 6.],
[5., 1., 4., 7.]])
# example with infinite eigenvalues
A4 = np.eye(2)
B4 = np.diag([0, 1])
# example with (alpha, beta) = (0, 0)
A5 = np.diag([1, 0])
cls.A = [A1, A2, A3, A4, A5]
cls.B = [B1, B2, B3, B4, A5]
def qz_decomp(self, sort):
with np.errstate(all='raise'):
ret = [ordqz(Ai, Bi, sort=sort) for Ai, Bi in zip(self.A, self.B)]
return tuple(ret)
def check(self, A, B, sort, AA, BB, alpha, beta, Q, Z):
Id = np.eye(*A.shape)
# make sure Q and Z are orthogonal
assert_array_almost_equal(Q @ Q.T.conj(), Id)
assert_array_almost_equal(Z @ Z.T.conj(), Id)
# check factorization
assert_array_almost_equal(Q @ AA, A @ Z)
assert_array_almost_equal(Q @ BB, B @ Z)
# check shape of AA and BB
assert_array_equal(np.tril(AA, -2), np.zeros(AA.shape))
assert_array_equal(np.tril(BB, -1), np.zeros(BB.shape))
# check eigenvalues
for i in range(A.shape[0]):
# does the current diagonal element belong to a 2-by-2 block
# that was already checked?
if i > 0 and A[i, i - 1] != 0:
continue
# take care of 2-by-2 blocks
if i < AA.shape[0] - 1 and AA[i + 1, i] != 0:
evals, _ = eig(AA[i:i + 2, i:i + 2], BB[i:i + 2, i:i + 2])
# make sure the pair of complex conjugate eigenvalues
# is ordered consistently (positive imaginary part first)
if evals[0].imag < 0:
evals = evals[[1, 0]]
tmp = alpha[i:i + 2]/beta[i:i + 2]
if tmp[0].imag < 0:
tmp = tmp[[1, 0]]
assert_array_almost_equal(evals, tmp)
else:
if alpha[i] == 0 and beta[i] == 0:
assert_equal(AA[i, i], 0)
assert_equal(BB[i, i], 0)
elif beta[i] == 0:
assert_equal(BB[i, i], 0)
else:
assert_almost_equal(AA[i, i]/BB[i, i], alpha[i]/beta[i])
sortfun = _select_function(sort)
lastsort = True
for i in range(A.shape[0]):
cursort = sortfun(np.array([alpha[i]]), np.array([beta[i]]))
# once the sorting criterion was not matched all subsequent
# eigenvalues also shouldn't match
if not lastsort:
assert not cursort
lastsort = cursort
def check_all(self, sort):
ret = self.qz_decomp(sort)
for reti, Ai, Bi in zip(ret, self.A, self.B):
self.check(Ai, Bi, sort, *reti)
def test_lhp(self):
self.check_all('lhp')
def test_rhp(self):
self.check_all('rhp')
def test_iuc(self):
self.check_all('iuc')
def test_ouc(self):
self.check_all('ouc')
def test_ref(self):
# real eigenvalues first (top-left corner)
def sort(x, y):
out = np.empty_like(x, dtype=bool)
nonzero = (y != 0)
out[~nonzero] = False
out[nonzero] = (x[nonzero]/y[nonzero]).imag == 0
return out
self.check_all(sort)
def test_cef(self):
# complex eigenvalues first (top-left corner)
def sort(x, y):
out = np.empty_like(x, dtype=bool)
nonzero = (y != 0)
out[~nonzero] = False
out[nonzero] = (x[nonzero]/y[nonzero]).imag != 0
return out
self.check_all(sort)
def test_diff_input_types(self):
ret = ordqz(self.A[1], self.B[2], sort='lhp')
self.check(self.A[1], self.B[2], 'lhp', *ret)
ret = ordqz(self.B[2], self.A[1], sort='lhp')
self.check(self.B[2], self.A[1], 'lhp', *ret)
def test_sort_explicit(self):
# Test order of the eigenvalues in the 2 x 2 case where we can
# explicitly compute the solution
A1 = np.eye(2)
B1 = np.diag([-2, 0.5])
expected1 = [('lhp', [-0.5, 2]),
('rhp', [2, -0.5]),
('iuc', [-0.5, 2]),
('ouc', [2, -0.5])]
A2 = np.eye(2)
B2 = np.diag([-2 + 1j, 0.5 + 0.5j])
expected2 = [('lhp', [1/(-2 + 1j), 1/(0.5 + 0.5j)]),
('rhp', [1/(0.5 + 0.5j), 1/(-2 + 1j)]),
('iuc', [1/(-2 + 1j), 1/(0.5 + 0.5j)]),
('ouc', [1/(0.5 + 0.5j), 1/(-2 + 1j)])]
# 'lhp' is ambiguous so don't test it
A3 = np.eye(2)
B3 = np.diag([2, 0])
expected3 = [('rhp', [0.5, np.inf]),
('iuc', [0.5, np.inf]),
('ouc', [np.inf, 0.5])]
# 'rhp' is ambiguous so don't test it
A4 = np.eye(2)
B4 = np.diag([-2, 0])
expected4 = [('lhp', [-0.5, np.inf]),
('iuc', [-0.5, np.inf]),
('ouc', [np.inf, -0.5])]
A5 = np.diag([0, 1])
B5 = np.diag([0, 0.5])
# 'lhp' and 'iuc' are ambiguous so don't test them
expected5 = [('rhp', [2, np.nan]),
('ouc', [2, np.nan])]
A = [A1, A2, A3, A4, A5]
B = [B1, B2, B3, B4, B5]
expected = [expected1, expected2, expected3, expected4, expected5]
for Ai, Bi, expectedi in zip(A, B, expected):
for sortstr, expected_eigvals in expectedi:
_, _, alpha, beta, _, _ = ordqz(Ai, Bi, sort=sortstr)
azero = (alpha == 0)
bzero = (beta == 0)
x = np.empty_like(alpha)
x[azero & bzero] = np.nan
x[~azero & bzero] = np.inf
x[~bzero] = alpha[~bzero]/beta[~bzero]
assert_allclose(expected_eigvals, x)
class TestOrdQZWorkspaceSize:
def test_decompose(self):
rng = np.random.RandomState(12345)
N = 202
# raises error if lwork parameter to dtrsen is too small
for ddtype in [np.float32, np.float64]:
A = rng.random((N, N)).astype(ddtype)
B = rng.random((N, N)).astype(ddtype)
# sort = lambda ar, ai, b: ar**2 + ai**2 < b**2
_ = ordqz(A, B, sort=lambda alpha, beta: alpha < beta,
output='real')
for ddtype in [np.complex128, np.complex64]:
A = rng.random((N, N)).astype(ddtype)
B = rng.random((N, N)).astype(ddtype)
_ = ordqz(A, B, sort=lambda alpha, beta: alpha < beta,
output='complex')
@pytest.mark.slow
def test_decompose_ouc(self):
rng = np.random.RandomState(12345)
N = 202
# segfaults if lwork parameter to dtrsen is too small
for ddtype in [np.float32, np.float64, np.complex128, np.complex64]:
A = rng.random((N, N)).astype(ddtype)
B = rng.random((N, N)).astype(ddtype)
S, T, alpha, beta, U, V = ordqz(A, B, sort='ouc')
class TestDatacopied:
def test_datacopied(self):
from scipy.linalg._decomp import _datacopied
M = matrix([[0, 1], [2, 3]])
A = asarray(M)
L = M.tolist()
M2 = M.copy()
class Fake1:
def __array__(self, dtype=None, copy=None):
return A
class Fake2:
__array_interface__ = A.__array_interface__
F1 = Fake1()
F2 = Fake2()
for item, status in [(M, False), (A, False), (L, True),
(M2, False), (F1, False), (F2, False)]:
arr = asarray(item)
assert_equal(_datacopied(arr, item), status,
err_msg=repr(item))
def test_aligned_mem_float():
"""Check linalg works with non-aligned memory (float32)"""
# Allocate 402 bytes of memory (allocated on boundary)
a = arange(402, dtype=np.uint8)
# Create an array with boundary offset 4
z = np.frombuffer(a.data, offset=2, count=100, dtype=float32)
z.shape = 10, 10
eig(z, overwrite_a=True)
eig(z.T, overwrite_a=True)
@pytest.mark.skipif(platform.machine() == 'ppc64le',
reason="crashes on ppc64le")
def test_aligned_mem():
"""Check linalg works with non-aligned memory (float64)"""
# Allocate 804 bytes of memory (allocated on boundary)
a = arange(804, dtype=np.uint8)
# Create an array with boundary offset 4
z = np.frombuffer(a.data, offset=4, count=100, dtype=float)
z.shape = 10, 10
eig(z, overwrite_a=True)
eig(z.T, overwrite_a=True)
def test_aligned_mem_complex():
"""Check that complex objects don't need to be completely aligned"""
# Allocate 1608 bytes of memory (allocated on boundary)
a = zeros(1608, dtype=np.uint8)
# Create an array with boundary offset 8
z = np.frombuffer(a.data, offset=8, count=100, dtype=complex)
z.shape = 10, 10
eig(z, overwrite_a=True)
# This does not need special handling
eig(z.T, overwrite_a=True)
def check_lapack_misaligned(func, args, kwargs):
args = list(args)
for i in range(len(args)):
a = args[:]
if isinstance(a[i], np.ndarray):
# Try misaligning a[i]
aa = np.zeros(a[i].size*a[i].dtype.itemsize+8, dtype=np.uint8)
aa = np.frombuffer(aa.data, offset=4, count=a[i].size,
dtype=a[i].dtype)
aa.shape = a[i].shape
aa[...] = a[i]
a[i] = aa
func(*a, **kwargs)
if len(a[i].shape) > 1:
a[i] = a[i].T
func(*a, **kwargs)
@pytest.mark.xfail(run=False,
reason="Ticket #1152, triggers a segfault in rare cases.")
def test_lapack_misaligned():
M = np.eye(10, dtype=float)
R = np.arange(100)
R.shape = 10, 10
S = np.arange(20000, dtype=np.uint8)
S = np.frombuffer(S.data, offset=4, count=100, dtype=float)
S.shape = 10, 10
b = np.ones(10)
LU, piv = lu_factor(S)
for (func, args, kwargs) in [
(eig, (S,), dict(overwrite_a=True)), # crash
(eigvals, (S,), dict(overwrite_a=True)), # no crash
(lu, (S,), dict(overwrite_a=True)), # no crash
(lu_factor, (S,), dict(overwrite_a=True)), # no crash
(lu_solve, ((LU, piv), b), dict(overwrite_b=True)),
(solve, (S, b), dict(overwrite_a=True, overwrite_b=True)),
(svd, (M,), dict(overwrite_a=True)), # no crash
(svd, (R,), dict(overwrite_a=True)), # no crash
(svd, (S,), dict(overwrite_a=True)), # crash
(svdvals, (S,), dict()), # no crash
(svdvals, (S,), dict(overwrite_a=True)), # crash
(cholesky, (M,), dict(overwrite_a=True)), # no crash
(qr, (S,), dict(overwrite_a=True)), # crash
(rq, (S,), dict(overwrite_a=True)), # crash
(hessenberg, (S,), dict(overwrite_a=True)), # crash
(schur, (S,), dict(overwrite_a=True)), # crash
]:
check_lapack_misaligned(func, args, kwargs)
# not properly tested
# cholesky, rsf2csf, lu_solve, solve, eig_banded, eigvals_banded, eigh, diagsvd
class TestOverwrite:
def test_eig(self):
assert_no_overwrite(eig, [(3, 3)])
assert_no_overwrite(eig, [(3, 3), (3, 3)])
def test_eigh(self):
assert_no_overwrite(eigh, [(3, 3)])
assert_no_overwrite(eigh, [(3, 3), (3, 3)])
def test_eig_banded(self):
assert_no_overwrite(eig_banded, [(3, 2)])
def test_eigvals(self):
assert_no_overwrite(eigvals, [(3, 3)])
def test_eigvalsh(self):
assert_no_overwrite(eigvalsh, [(3, 3)])
def test_eigvals_banded(self):
assert_no_overwrite(eigvals_banded, [(3, 2)])
def test_hessenberg(self):
assert_no_overwrite(hessenberg, [(3, 3)])
def test_lu_factor(self):
assert_no_overwrite(lu_factor, [(3, 3)])
def test_lu_solve(self):
x = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 8]])
xlu = lu_factor(x)
assert_no_overwrite(lambda b: lu_solve(xlu, b), [(3,)])
def test_lu(self):
assert_no_overwrite(lu, [(3, 3)])
def test_qr(self):
assert_no_overwrite(qr, [(3, 3)])
def test_rq(self):
assert_no_overwrite(rq, [(3, 3)])
def test_schur(self):
assert_no_overwrite(schur, [(3, 3)])
def test_schur_complex(self):
assert_no_overwrite(lambda a: schur(a, 'complex'), [(3, 3)],
dtypes=[np.float32, np.float64])
def test_svd(self):
assert_no_overwrite(svd, [(3, 3)])
assert_no_overwrite(lambda a: svd(a, lapack_driver='gesvd'), [(3, 3)])
def test_svdvals(self):
assert_no_overwrite(svdvals, [(3, 3)])
def _check_orth(n, dtype, skip_big=False):
X = np.ones((n, 2), dtype=float).astype(dtype)
eps = np.finfo(dtype).eps
tol = 1000 * eps
Y = orth(X)
assert_equal(Y.shape, (n, 1))
assert_allclose(Y, Y.mean(), atol=tol)
Y = orth(X.T)
assert_equal(Y.shape, (2, 1))
assert_allclose(Y, Y.mean(), atol=tol)
if n > 5 and not skip_big:
np.random.seed(1)
X = np.random.rand(n, 5) @ np.random.rand(5, n)
X = X + 1e-4 * np.random.rand(n, 1) @ np.random.rand(1, n)
X = X.astype(dtype)
Y = orth(X, rcond=1e-3)
assert_equal(Y.shape, (n, 5))
Y = orth(X, rcond=1e-6)
assert_equal(Y.shape, (n, 5 + 1))
@pytest.mark.slow
@pytest.mark.skipif(np.dtype(np.intp).itemsize < 8,
reason="test only on 64-bit, else too slow")
def test_orth_memory_efficiency():
# Pick n so that 16*n bytes is reasonable but 8*n*n bytes is unreasonable.
# Keep in mind that @pytest.mark.slow tests are likely to be running
# under configurations that support 4Gb+ memory for tests related to
# 32 bit overflow.
n = 10*1000*1000
try:
_check_orth(n, np.float64, skip_big=True)
except MemoryError as e:
raise AssertionError(
'memory error perhaps caused by orth regression'
) from e
def test_orth():
dtypes = [np.float32, np.float64, np.complex64, np.complex128]
sizes = [1, 2, 3, 10, 100]
for dt, n in itertools.product(dtypes, sizes):
_check_orth(n, dt)
def test_null_space():
np.random.seed(1)
dtypes = [np.float32, np.float64, np.complex64, np.complex128]
sizes = [1, 2, 3, 10, 100]
for dt, n in itertools.product(dtypes, sizes):
X = np.ones((2, n), dtype=dt)
eps = np.finfo(dt).eps
tol = 1000 * eps
Y = null_space(X)
assert_equal(Y.shape, (n, n-1))
assert_allclose(X @ Y, 0, atol=tol)
Y = null_space(X.T)
assert_equal(Y.shape, (2, 1))
assert_allclose(X.T @ Y, 0, atol=tol)
X = np.random.randn(1 + n//2, n)
Y = null_space(X)
assert_equal(Y.shape, (n, n - 1 - n//2))
assert_allclose(X @ Y, 0, atol=tol)
if n > 5:
np.random.seed(1)
X = np.random.rand(n, 5) @ np.random.rand(5, n)
X = X + 1e-4 * np.random.rand(n, 1) @ np.random.rand(1, n)
X = X.astype(dt)
Y = null_space(X, rcond=1e-3)
assert_equal(Y.shape, (n, n - 5))
Y = null_space(X, rcond=1e-6)
assert_equal(Y.shape, (n, n - 6))
def test_subspace_angles():
H = hadamard(8, float)
A = H[:, :3]
B = H[:, 3:]
assert_allclose(subspace_angles(A, B), [np.pi / 2.] * 3, atol=1e-14)
assert_allclose(subspace_angles(B, A), [np.pi / 2.] * 3, atol=1e-14)
for x in (A, B):
assert_allclose(subspace_angles(x, x), np.zeros(x.shape[1]),
atol=1e-14)
# From MATLAB function "subspace", which effectively only returns the
# last value that we calculate
x = np.array(
[[0.537667139546100, 0.318765239858981, 3.578396939725760, 0.725404224946106], # noqa: E501
[1.833885014595086, -1.307688296305273, 2.769437029884877, -0.063054873189656], # noqa: E501
[-2.258846861003648, -0.433592022305684, -1.349886940156521, 0.714742903826096], # noqa: E501
[0.862173320368121, 0.342624466538650, 3.034923466331855, -0.204966058299775]]) # noqa: E501
expected = 1.481454682101605
assert_allclose(subspace_angles(x[:, :2], x[:, 2:])[0], expected,
rtol=1e-12)
assert_allclose(subspace_angles(x[:, 2:], x[:, :2])[0], expected,
rtol=1e-12)
expected = 0.746361174247302
assert_allclose(subspace_angles(x[:, :2], x[:, [2]]), expected, rtol=1e-12)
assert_allclose(subspace_angles(x[:, [2]], x[:, :2]), expected, rtol=1e-12)
expected = 0.487163718534313
assert_allclose(subspace_angles(x[:, :3], x[:, [3]]), expected, rtol=1e-12)
assert_allclose(subspace_angles(x[:, [3]], x[:, :3]), expected, rtol=1e-12)
expected = 0.328950515907756
assert_allclose(subspace_angles(x[:, :2], x[:, 1:]), [expected, 0],
atol=1e-12)
# Degenerate conditions
assert_raises(ValueError, subspace_angles, x[0], x)
assert_raises(ValueError, subspace_angles, x, x[0])
assert_raises(ValueError, subspace_angles, x[:-1], x)
# Test branch if mask.any is True:
A = np.array([[1, 0, 0],
[0, 1, 0],
[0, 0, 1],
[0, 0, 0],
[0, 0, 0]])
B = np.array([[1, 0, 0],
[0, 1, 0],
[0, 0, 0],
[0, 0, 0],
[0, 0, 1]])
expected = np.array([np.pi/2, 0, 0])
assert_allclose(subspace_angles(A, B), expected, rtol=1e-12)
# Complex
# second column in "b" does not affect result, just there so that
# b can have more cols than a, and vice-versa (both conditional code paths)
a = [[1 + 1j], [0]]
b = [[1 - 1j, 0], [0, 1]]
assert_allclose(subspace_angles(a, b), 0., atol=1e-14)
assert_allclose(subspace_angles(b, a), 0., atol=1e-14)
class TestCDF2RDF:
def matmul(self, a, b):
return np.einsum('...ij,...jk->...ik', a, b)
def assert_eig_valid(self, w, v, x):
assert_array_almost_equal(
self.matmul(v, w),
self.matmul(x, v)
)
def test_single_array0x0real(self):
# eig doesn't support 0x0 in old versions of numpy
X = np.empty((0, 0))
w, v = np.empty(0), np.empty((0, 0))
wr, vr = cdf2rdf(w, v)
self.assert_eig_valid(wr, vr, X)
def test_single_array2x2_real(self):
X = np.array([[1, 2], [3, -1]])
w, v = np.linalg.eig(X)
wr, vr = cdf2rdf(w, v)
self.assert_eig_valid(wr, vr, X)
def test_single_array2x2_complex(self):
X = np.array([[1, 2], [-2, 1]])
w, v = np.linalg.eig(X)
wr, vr = cdf2rdf(w, v)
self.assert_eig_valid(wr, vr, X)
def test_single_array3x3_real(self):
X = np.array([[1, 2, 3], [1, 2, 3], [2, 5, 6]])
w, v = np.linalg.eig(X)
wr, vr = cdf2rdf(w, v)
self.assert_eig_valid(wr, vr, X)
def test_single_array3x3_complex(self):
X = np.array([[1, 2, 3], [0, 4, 5], [0, -5, 4]])
w, v = np.linalg.eig(X)
wr, vr = cdf2rdf(w, v)
self.assert_eig_valid(wr, vr, X)
def test_random_1d_stacked_arrays(self):
# cannot test M == 0 due to bug in old numpy
for M in range(1, 7):
np.random.seed(999999999)
X = np.random.rand(100, M, M)
w, v = np.linalg.eig(X)
wr, vr = cdf2rdf(w, v)
self.assert_eig_valid(wr, vr, X)
def test_random_2d_stacked_arrays(self):
# cannot test M == 0 due to bug in old numpy
for M in range(1, 7):
X = np.random.rand(10, 10, M, M)
w, v = np.linalg.eig(X)
wr, vr = cdf2rdf(w, v)
self.assert_eig_valid(wr, vr, X)
def test_low_dimensionality_error(self):
w, v = np.empty(()), np.array((2,))
assert_raises(ValueError, cdf2rdf, w, v)
def test_not_square_error(self):
# Check that passing a non-square array raises a ValueError.
w, v = np.arange(3), np.arange(6).reshape(3, 2)
assert_raises(ValueError, cdf2rdf, w, v)
def test_swapped_v_w_error(self):
# Check that exchanging places of w and v raises ValueError.
X = np.array([[1, 2, 3], [0, 4, 5], [0, -5, 4]])
w, v = np.linalg.eig(X)
assert_raises(ValueError, cdf2rdf, v, w)
def test_non_associated_error(self):
# Check that passing non-associated eigenvectors raises a ValueError.
w, v = np.arange(3), np.arange(16).reshape(4, 4)
assert_raises(ValueError, cdf2rdf, w, v)
def test_not_conjugate_pairs(self):
# Check that passing non-conjugate pairs raises a ValueError.
X = np.array([[1, 2, 3], [1, 2, 3], [2, 5, 6+1j]])
w, v = np.linalg.eig(X)
assert_raises(ValueError, cdf2rdf, w, v)
# different arrays in the stack, so not conjugate
X = np.array([
[[1, 2, 3], [1, 2, 3], [2, 5, 6+1j]],
[[1, 2, 3], [1, 2, 3], [2, 5, 6-1j]],
])
w, v = np.linalg.eig(X)
assert_raises(ValueError, cdf2rdf, w, v)