From 8e621e14a86d457dcae6f78f6b37263995382d51 Mon Sep 17 00:00:00 2001 From: Tim Widrick Date: Sun, 1 Oct 2023 11:06:05 -0400 Subject: [PATCH] update ytools.mattype: add option to output cholesky decomp and better type checking update ytools.eig_si to request cholesky decomp from mattype (faster) --- pyyeti/tests/test_ytools.py | 32 +++++++++- pyyeti/ytools.py | 116 ++++++++++++++++++++++-------------- 2 files changed, 103 insertions(+), 45 deletions(-) diff --git a/pyyeti/tests/test_ytools.py b/pyyeti/tests/test_ytools.py index 7c1a4095..0d179b3c 100644 --- a/pyyeti/tests/test_ytools.py +++ b/pyyeti/tests/test_ytools.py @@ -3,7 +3,7 @@ import tempfile import os from pyyeti import ytools -import scipy.linalg as linalg +from scipy import linalg import pytest @@ -130,6 +130,36 @@ def test_mattype(): ytools.mattype(c, "badtype") +def test_mattype2(): + # real a + a = np.array([[346500.0, 1e-7], [2.1e-7, 1000000.1]]) + assert 13 == ytools.mattype(a)[0] + t, mattypes, ch = ytools.mattype(a, return_cholesky=True) + assert t == 13 + assert np.allclose(ch, linalg.cholesky(a)) + + # complex a + a = np.array([[346500.0, 1e-7 * (1 + 1j)], [2.1e-7 * (1 - 1j), 1000000.1]]) + assert 15 == ytools.mattype(a)[0] + t, mattypes, ch = ytools.mattype(a, return_cholesky=True) + assert t == 15 + assert np.allclose(ch, linalg.cholesky(a)) + + # real a + a = np.array([[346500.0, 1e-7], [2000.1e-7, 1000000.1]]) + assert 0 == ytools.mattype(a)[0] + t, mattypes, ch = ytools.mattype(a, return_cholesky=True) + assert t == 0 + assert ch is None + + # complex a + a = np.array([[346500.0, 1e-7 * (1 + 1j)], [2000.1e-7 * (1 - 1j), 1000000.1]]) + assert 0 == ytools.mattype(a)[0] + t, mattypes, ch = ytools.mattype(a, return_cholesky=True) + assert t == 0 + assert ch is None + + def test_save_load(): a = np.arange(18).reshape(2, 3, 3) b = np.arange(3) diff --git a/pyyeti/ytools.py b/pyyeti/ytools.py index 46e7bc44..c0b79855 100644 --- a/pyyeti/ytools.py +++ b/pyyeti/ytools.py @@ -715,7 +715,40 @@ def isdiag(A, tol=1e-12): return max_off <= tol * max_on -def mattype(A, mtype=None): +def _check_symm_herm(A, mattypes): + Atype = 0 + if np.allclose(A, A.T): + Atype |= mattypes["symmetric"] + + elif np.iscomplexobj(A) and np.allclose(A, A.T.conj()): + Atype |= mattypes["hermitian"] + + if isdiag(A): + Atype |= mattypes["diagonal"] | mattypes["symmetric"] + if np.iscomplexobj(A): + Atype |= mattypes["hermitian"] + d = np.diag(A) + if np.allclose(1, d): + Atype |= mattypes["identity"] + + return Atype + + +def _check_cholesky(Atype, A, mattypes): + chol = None + if (Atype & mattypes["symmetric"] and np.isrealobj(A)) or ( + Atype & mattypes["hermitian"] + ): + try: + chol = linalg.cholesky(A) + except linalg.LinAlgError: + pass + else: + Atype |= mattypes["posdef"] + return Atype, chol + + +def mattype(A, mtype=None, return_cholesky=False): """ Checks contents of square matrix `A` to see if it is symmetric, hermitian, positive-definite, diagonal, and identity. @@ -731,6 +764,10 @@ def mattype(A, mtype=None): this case, True is returned if `A` is of the type specified or False otherwise. If None, `Atype` (if `A` is not None) and `mattypes` is returned. `mtype` is ignored if `A` is None. + return_cholesky : bool; optional + If True, the output of :func:`scipy.linalg.cholesky` is + returned if computed. Output will be None if `A` is not + positive-definite. See example usages below. Returns ------- @@ -753,13 +790,25 @@ def mattype(A, mtype=None): Not returned if `mtype` is specified. This is the only return if `A` is None. + chol : 2d ndarray or None + See `return_cholesky` above. If returned, `chol` will be the + output of :func:`scipy.linalg.cholesky` (with default + settings) or None, depending on whether matrix is + positive-definite or not. + Notes ----- - Here are some example usages:: - - mattype(A) # returns (Atype, mattypes) - mattype(A, 'symmetric') # returns True or False - mattype(None) # returns mattypes + Here are some example usages: + + ========================================== ======================= + Usage Returns + ========================================== ======================= + mattype(A) (Atype, mattypes) + mattype(A, return_cholesky=True) (Atype, mattypes, chol) + mattype(A, 'symmetric') True or False + mattype(A, 'posdef', return_cholesky=True) (True or False, chol) + mattype(None) mattypes + ========================================== ======================= See also -------- @@ -809,51 +858,30 @@ def mattype(A, mtype=None): if mtype is None: if A.ndim != 2 or A.shape[0] != A.shape[1]: return Atype, mattypes - if np.allclose(A, A.T): - Atype |= mattypes["symmetric"] - if np.isrealobj(A): - try: - linalg.cholesky(A) - except linalg.LinAlgError: - pass - else: - Atype |= mattypes["posdef"] - elif np.iscomplexobj(A) and np.allclose(A, A.T.conj()): - Atype |= mattypes["hermitian"] - try: - linalg.cholesky(A) - except linalg.LinAlgError: - pass - else: - Atype |= mattypes["posdef"] - if isdiag(A): - Atype |= mattypes["diagonal"] - d = np.diag(A) - if np.allclose(1, d): - Atype |= mattypes["identity"] + + Atype = _check_symm_herm(A, mattypes) + Atype, chol = _check_cholesky(Atype, A, mattypes) + + if return_cholesky: + return Atype, mattypes, chol return Atype, mattypes if A.ndim != 2 or A.shape[0] != A.shape[1]: return False if mtype == "symmetric": - return np.allclose(A, A.T) + return np.allclose(A, A.T) or isdiag(A) if mtype == "hermitian": - return np.allclose(A, A.T.conj()) + return np.iscomplexobj(A) and (np.allclose(A, A.T.conj()) or isdiag(A)) if mtype == "posdef": - if np.isrealobj(A): - if not np.allclose(A, A.T): - return False - else: - if not np.allclose(A, A.T.conj()): - return False - try: - linalg.cholesky(A) - return True - except linalg.LinAlgError: - return False + Atype = _check_symm_herm(A, mattypes) + Atype, chol = _check_cholesky(Atype, A, mattypes) + ret = bool(Atype & mattypes["posdef"]) + if return_cholesky: + return ret, chol + return ret if mtype in ("diagonal", "identity"): if isdiag(A): @@ -1148,7 +1176,7 @@ def eig_si( # Mk = (Mk + Mk.T) / 2 # solve subspace eigenvalue problem: - mtp = mattype(Mk, "posdef") + mtp, Mkuu = mattype(Mk, "posdef", return_cholesky=True) if not mtp: factor = 1000 * eps pc = 0 @@ -1156,12 +1184,12 @@ def eig_si( pc += 1 Mk += np.diag(np.diag(Mk) * factor) factor *= 10.0 - mtp = mattype(Mk, "posdef") + mtp, Mkuu = mattype(Mk, "posdef", return_cholesky=True) if mtp or pc > 5: break if mtp: - Mkll = linalg.cholesky(Mk, lower=False).T + Mkll = Mkuu.T # linalg.cholesky(Mk, lower=False).T Kkmod = linalg.solve_triangular( Mkll, linalg.solve_triangular(Mkll, Kk, lower=True).T, lower=True )