Source code for qsymm.linalg

import warnings
import numpy as np
import scipy.linalg as la
import scipy.sparse.linalg as sla
import scipy
from numbers import Number
from scipy.spatial.distance import cdist
from scipy.sparse.csgraph import connected_components
import itertools as it
import tinyarray as ta


[docs] def commutator(A, B): return np.dot(A, B) - np.dot(B, A)
[docs] def prop_to_id(A): # Test if A is proportional to the identity matrix # and return the factor as well if np.isclose(A[0, 0], 0): if np.allclose(A, np.zeros(A.shape)): return True, 0 else: return False, 0 else: if np.allclose(A / A[0, 0], np.eye(*A.shape)): return (True, A[0, 0]) else: return False, 0
[docs] def allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False): """Alternative to numpy.allclose to test that two ndarrays are elementwise close. Unlike the numpy version, the relative tolerance is not elementwise, but it is relative to the largest entry in the arrays.""" # If either is a list or number, recast it to ndarray, if it is # a LinearOperator, recast it to dense matrix. This is rather # inefficient, but LinearOperator doesn't support multiplication # with sparse matrices. if isinstance(a, (list, Number)): a = np.asarray(a) elif isinstance(a, sla.LinearOperator): a = a @ np.eye(a.shape[-1], dtype=complex) if isinstance(b, (list, Number)): b = np.asarray(b) elif isinstance(b, sla.LinearOperator): b = b @ np.eye(a.shape[-1], dtype=complex) # Check if empty arrays, only compare shape if a.size == 0: return a.shape == b.shape atol = atol + rtol * max(np.max(np.abs(a)), np.max(np.abs(b))) diff = a - b if isinstance(diff, scipy.sparse.spmatrix): diff = diff.data return np.allclose(diff, 0, rtol=0, atol=atol, equal_nan=equal_nan)
[docs] def mtm(a, B, c): # matrix-tensor-matrix multiplication for dimensions 2, 3, 2. # Equivalent to 'np.einsum('ij,ajk,kl->ail', a, B, c)' # and to 'a @ B @ c', but much faster. return np.array([a.dot(b).dot(c) for b in B])
[docs] def matrix_basis(dim, traceless=False, antihermitian=False, real=False, sparse=False): """ "Construct a basis for the vector space of dim x dim matrices, that may in addition be traceless, and/or composed of Hermitian or skew-Hermitian basis elements. Parameters ---------- dim: positive integer The dimension of the matrices (dim x dim). traceless: boolean Whether the vector space only includes traceless matrices (True) or not (False). antihermitian: boolean Whether to use a basis composed of Hermitian matrices (False), or skew-Hermitian matrices (True). real: Whether the basis only spans real matrices. Depending on 'antihermitian' it spans symmetric (False) or antisymmetric (True) matrices only. sparse: If True, scipy.sparse.csr_matrix is returned, otherwise np.ndarray. Returns --------- basis: generator A generator that returns matrices that span the vector space. """ if sparse: def null(): return scipy.sparse.lil_matrix((dim, dim), dtype=complex) def set_type(x): return x.tocsr() else: def null(): return np.zeros((dim, dim), dtype=complex) def set_type(x): return x # Matrix basis for dim x dim matrices. With real coefficients spans # Hermitian matrices, with complex spans all matrices coeff = 1j if antihermitian else 1 # Diagonals if real and antihermitian: pass elif traceless: for i in range(dim - 1): diag = null() diag[i, i] = 1 diag[dim - 1, dim - 1] = -1 yield coeff * set_type(diag) else: for i in range(dim): diag = null() diag[i, i] = 1 yield coeff * set_type(diag) # Off-diagonals for i, j in it.combinations(list(range(dim)), 2): if not (real and antihermitian): h = null() h[i, j] = 1 h[j, i] = 1 yield coeff * set_type(h) if not (real and not antihermitian): h = null() h[i, j] = 1j h[j, i] = -1j yield coeff * set_type(h)
[docs] def family_to_vectors(family, all_keys=None): """Convert a list of Model to standard vector representation. Parameters ---------- family: iterable of Model objects The Model to convert to vector representation. all_keys: iterable or None Iterable that must contain all the keys in the family. An optional parameter that allows to fix the ordering of the returned vectors to the same as in all_keys. Otherwise, the keys are ordered arbitrarily (all_keys = None). Returns ---------- A matrix representing the family, in which the family members form the rows.""" def matrix_to_vector(M): M = M.reshape(-1) return np.concatenate((M.real, M.imag)) if all_keys is None: all_keys = set() for term in family: all_keys |= term.keys() all_keys = list(all_keys) assert all([set(all_keys) >= term.keys() for term in family]) if not all_keys: return np.empty((len(family), 0)) # Map each family member to a vector, gather the vectors into # a matrix in which the vectors form the rows. basis_vectors = [] for member in family: vector = [matrix_to_vector(mat) for mat in member.value_list(all_keys)] basis_vectors.append(np.hstack(vector)) return np.vstack(basis_vectors)
[docs] def nullspace(A, atol=1e-6, return_complement=False, sparse=None, k_max=-10): """Compute an approximate basis for the nullspace of matrix A. Parameters ---------- A : ndarray or spmatrix A should be 2-D. atol : real Absolute tolerance when deciding whether an eigen/singular value is zero, so the corresponding vector is included in the null space. return_complement : bool Whether to return the basis of the orthocomplement ot the null space as well. Treated as False if sparse is True. sparse : bool or None (default) Whether to use sparse linear algebra to find the null space. The default value is inferred from the type of the input array. k_max : int (default -10) If k_max<0, the full null-space is found by increasing the number of eigenvectors found in each step by abs(k_max). If k_max>0, k_max is the maximum number of null space vectors requested. If the null space dimension is higher than k_max, basis of a random subspace is returned. Can be used to increase performance if the maximum null space dimensionality is known. Ignored unless sparse is True. Returns ------- ns : ndarray Basis of the null space of A as row vectors. nsc : ndarray Basis of the complement of the null space as row vectors, the othonormal basis of the row span of A. Only returned if return_complement=True and sparse=False. """ if sparse is None: sparse = isinstance(A, scipy.sparse.spmatrix) if sparse: # Do sparse eigenvalue solving # sigma used for shift-invert, should be a small negative value sigma = -1e-1 # Make sure A is sparse matrix if isinstance(A, np.ndarray): A = scipy.sparse.csr_matrix(A) A.eliminate_zeros() # Treat A=0 case if np.allclose(A.data, 0, atol=atol / A.shape[1]): return np.eye(A.shape[1]) # Make Hermitian positive definite matrix A = A.T.conj().dot(A) if k_max > 0: # If k_max is specified, find k_max smallest eigenvalues evals, evecs = sla.eigsh( A, sigma=sigma, which="LM", k=k_max, return_eigenvectors=True ) elif k_max < 0: # Successively find more eigenvalues until some not close to zero # Number of new null space vectors to find in each step k_step = abs(k_max) k_max = min(A.shape[0] - 2, k_step) while k_max < A.shape[0] - 1: evals, evecs = sla.eigsh( A, sigma=sigma, which="LM", k=k_max, return_eigenvectors=True ) if np.max(np.abs(evals)) > atol: # We found the first large one break else: # All small, need to increase k_max k_max += k_step else: raise ValueError( "A should have at most A.shape[0]-2 dimensional null space." "Try using sparse=False." ) else: raise ValueError("k_max must be nonzero.") # Only keep eigenvectors with small eigenvalues nnz = np.isclose(evals, 0, atol=atol) ns = evecs[:, nnz] # Orthonormalize null space if ns.shape[1] > 0: ns, _ = la.qr(ns, mode="economic") return ns.T else: if scipy.sparse.issparse(A): A = A.toarray() _, s, vh = la.svd(A, full_matrices=True) # Count singular values <= atol. s is sorted descending; reverse to get ascending. large_count = len(s) - np.searchsorted(s[::-1], atol, side="right") ns = vh[large_count:].conj() if return_complement: nsc = vh[:large_count].conj() return ns, nsc else: return ns
[docs] def split_list(vals, tol=1e-6): # Returns start and end indices of blocks of # consecutive close values in a list boundaries = np.where(np.abs(np.diff(vals)) > tol / len(vals))[0] + 1 return np.array([np.insert(boundaries, 0, 0), np.append(boundaries, len(vals))]).T
[docs] def grouped_diag(H, tol=1e-6, checks=0): r""" Diagonalize normal matrix 'H' and group the eigenvalues and eigenvectors such that approximately equal eigenvalues are grouped together. If 'H' is Hermitian, 'eigh' returns ordered eigenvalues. Returns the grouped eigenvalues, eigenvectors. """ Hdag = H.T.conj() if allclose(H, Hdag): evals, U = la.eigh(H) return evals, U if not np.allclose(commutator(H, Hdag), 0): raise ValueError("Only normal matrix can be diagonalized.") evals, U = la.eig(H) # `eig` does not order nearly degenerate complex eigenvalues reliably, so # cluster them in the complex plane before recursing on the blocks. evvec = np.array([evals.real, evals.imag]).T # Find connected clusters of close values con = cdist(evvec, evvec) < tol / len(H) _, groups = connected_components(con) # reorder evals and evecs such that groups are together order = np.argsort(groups) evals = evals[order] U = U[:, order] # QR only rotates vectors inside each clustered eigenspace, which is enough # to recover a unitary basis for the later block recursion. U, _ = la.qr(U) if checks == 2: # Check the result assert allclose((U.T.conj() @ H @ U), np.diag(evals)) return evals, U
[docs] def simult_diag(mats, tol=1e-6, checks=0): """ Simultaneously diagonalize commuting normal square matrices. Recursive algorithm, first diagonalize first matrix, express the rest of the matrices in its eigenbasis, and repeat the procedure for the rest of the matrices in each of the degenerate eigensubspaces. Recursion depth is limited by the number of nontrivial nested eigensubspaces, which is at most the number of matrices or the size of the matrices. Parameters ---------- mats : ndarray List of commuting matrices, 3 index array. tol : float Tolerance when deciding degeneracy of eigenvalues. checks : int Amount of checks to do. 0: no checks 1: check the final result 2: check all the intermediate results Returns ------- Ps : list of ndarrays List of rectangular matrices of common eigenvectors spanning each of the common eigensubspaces. Stacking them gives unitary matrix 'U = np.hstack(Ps)' such that 'U.T.conj() @ mats[i] @ U' is diagonal for every 'i'. Notes: It is recommended to order matrices such that those with eigenvalues that are far spaced or exactly degenerate (such as symmetries) are first, and those which may have accidental near degeneracies (such as Hamiltonians) are last. Tolerance is divided by the size of the block when deciding which eigenvalues are degenerate. The default value seems to work well, but not extensively tested, numerical instabilities are possible. """ # If 1x1 matrix, we are done if mats[0].shape == (1, 1): return [np.eye(1)] # Diagonalize mats[0], if there are no more matrices we are done evals, U = grouped_diag(mats[0], tol=tol, checks=checks) ind = split_list(evals, tol=tol) if len(mats) == 1: return [U[:, b:e] for b, e in ind] # Check that all matrices commute with mats[0] if not np.allclose([commutator(mats[0], mat) for mat in mats[1:]], 0): raise ValueError("Only commuting matrices can be simultaneously diagonalized.") # Transform the rest of the matrices mats[1:] to the eigenbasis # of mats[0] where they are all block-diagonal. # Apply the same algorithm to the rest of the matrices in # each block recursively. matr = mtm(U.T.conj(), mats[1:], U) if checks == 2: # Check that off-diagonal blocks are small after the transformation assert np.all( [ np.allclose(matr[:, b1:e1, b2:e2], 0) for (b1, e1), (b2, e2) in it.combinations(ind, 2) ] ) Ps = [] for b, e in ind: P0 = U[:, b:e] Pnew = simult_diag(matr[:, b:e, b:e], tol=tol, checks=(2 if checks == 2 else 0)) Ps += [np.dot(P0, P) for P in Pnew] if checks > 0: # Check the result is diagonal U = np.hstack(Ps) matsd = mtm(U.T.conj(), mats, U) assert np.all([allclose(m, np.diag(np.diagonal(m))) for m in matsd]) return Ps
[docs] def symmetry_adapted_sun(gens, check=False, n=None): r""" Find symmetry adapted basis of the simple 'su(d)' Lie-algebra representation defined by generators 'gens'. It is assumed that the representation is the direct sum of 'n' identical 'su(d)' representations, i.e. in the symmetry adapted basis the generators have the form 'L \otimes 1_{nxn}' where 'L' runs over a generator set of all 'd*d' Hermitian matrices. Also works for a direct sum of 'n' identical irreducible representations, in this case 'n' must be provided. Parameters ---------- gens : ndarray 3 index array of shape '(k, n*d, n*d)' with 'n' and 'd' integers, list of group generators. If 'n' is not provided, it is assumed that 'k = d**2-1' and it is a list of 'su(d)' generators. n : int or None (default) Number of copies of the identical irreducible representation. check : bool Whether to check the final result. Returns ------- Ps : ndarray 3 index array of shape '(d, n*d, n)', list of projectors onto the symmetry adapted subspaces. 'np.einsum('aij,ab,bkj->ik', Ps, L, Ps.conjugate()))' is a symmetry spanned by 'gens' for any 'd*d' Hermitian matrix 'L'. Stacking these produces the unitary transformation 'U = np.hstack(Ps)' to the symmetry adapted basis. """ if n is None: d = np.sqrt(len(gens) + 1) n = gens.shape[-1] / d if not (n.is_integer() and d.is_integer()): raise ValueError( "Shape of gens is incompatible with it being direct sum of" "n identical copies of su(d) representation." ) n, d = int(n), int(d) # Trivial case when it is a full SU(N) ### TODO: maybe the relative phase fixing would be still useful, # tests pass with this removed if n == 1: return np.array(np.split(np.eye(n * d), n * d, 1)) else: d = gens.shape[2] / n if not d.is_integer(): raise ValueError( "Shape of gens is incompatible with it being direct sum of" "n identical copies of an irreducible representation." ) d = int(d) # Iteratively split the irreps. Goal is to find a basis in the n*d dimensional # Hilbert-space, where the n vectors belonging to the identity factor in the symmetry # adapted basis are grouped together. A degenerate eigensubspace of a generator has # dimension of a multiple of n. If it is exactly n, we can be sure that we found # such a subspace. If it is f*n, it is the union of f such subspaces. We project the # remaining generators in this subspace, they will still span all matrices with the # structure L \otimes 1_{nxn} with L f*f. We pick a new generator and repeat the process. # In the generic case (a random combination of the generators) the first generator # already splits all the irreps, but in the standard matrix basis this is not the case. unsplit = [np.eye(n * d)] split = [] for g in gens: Ps = unsplit unsplit = [] for P in Ps: # project generator in unsplit subspace, rest of the generators still # span the full space of Hermitian matrices in the restricted space gr = P.T.conj() @ g @ P evals, U = grouped_diag(gr) assert allclose(U @ U.T.conj(), np.eye(U.shape[0])), U # find the degenerate eigenspaces blocks = split_list(evals) for b, e in blocks: if e - b == n: # if n long degenerate block, it is split split.append(P @ U[:, b:e]) else: assert (e - b) % n == 0 # if it is a multiple of n, put it back in unsplit unsplit.append(P @ U[:, b:e]) if unsplit == []: break else: # we could not split all the subspaces raise ValueError( "Algorithm failed, likely not a direct sum of identical irreducible representations." ) Ps = np.array(split) # Iterative phase fixing, start with first block fixed but unused fixed_unused = [0] fixed_used = [] unfixed = list(range(1, len(Ps))) while unfixed: # Adjust the basis of every block to the basis of a fixed block, # such that all the generators are proportional to the identity # in every block. try: j = fixed_unused.pop() except IndexError: # we could not fix some of the bases raise ValueError( "Algorithm failed, likely not a direct sum of identical irreducible representations." ) fixed_used.append(j) P0 = Ps[j] for i in unfixed: P = Ps[i] for g in gens: # find a generator where the block is nonzero gblock = P0.T.conj() @ g @ P if np.allclose(gblock, 0): continue # then it is proportional to unitary, transform it # to proportional to identity u, s, vh = la.svd(gblock) U = ((u @ vh).T).conj() Ps[i] = Ps[i] @ U unfixed.remove(i) fixed_unused.append(i) break if check: # check that every block is proportional to the identity U = np.hstack(Ps) genst = mtm(U.T.conjugate(), gens, U) assert np.all( [ [ prop_to_id(genst[i, n * j1 : n * (j1 + 1), n * j2 : n * (j2 + 1)])[ 0 ] for j1, j2 in it.combinations_with_replacement(range(d), 2) ] for i in range(len(genst)) ] ) return Ps
[docs] def solve_mat_eqn( HL, HR=None, hermitian=False, traceless=False, conjugate=False, sparse=None, k_max=-10, ): """Solve for X the simultaneous matrix equatioins X HL[i] = HR[i] X for every i. It is mapped to a system of linear equations, the null space of which gives a basis for all sulutions. Parameters ---------- HL : ndarray or list of arrays or sparse matrices Coefficient matrices of identical square shape, list of arrays of shape (n, n) or one array of shape (m, n, n). HR : ndarray or list of ndarrays or None Same as HL. If HR = None, HR = HL is used. hermitian, traceless : booleans If X has to be hermitian, use hermitian = True. If X has to be traceless, use traceless = True. conjugate : boolean or list of booleans If True, solve X HL[i] = HR[i] X^* instead. If a list with the same length as HL and HR is provided, conjugation is applied to the equations with index corresponding to the True entries. sparse : bool or None (default) Whether to use sparse linear algebra to find the null space. The default value is inferred from the type of the input array. k_max : int If k_max<0, all solutions are found by increasing the number of soulutions sought in each step by abs(k_max). If k_max>0, k_max is the maximum number of solutions requested. If the solution space dimension is higher than k_max, basis of a random subspace is returned. Can be used to increase performance if the maximum number of solutions is known. Ignored unless sparse is True. Returns ------- ndarray of shape (l, n, n), or list of csr_matrix List of linearly independent square matrices that span all solutions of the eaquations. """ if sparse is None: sparse = isinstance(HL[0], scipy.sparse.spmatrix) if HR is None: HR = HL if len(HL) != len(HR) or not all( left.shape == right.shape for left, right in zip(HL, HR) ): raise ValueError("HL and HR must have the same shape.") if isinstance(conjugate, bool): conjugate = [conjugate] * len(HL) if len(conjugate) != len(HL): raise ValueError("Conjugate must have the same length as HL.") if not all(term.shape[0] == term.shape[1] for term in HL): raise ValueError("HL and HR must be a list of square matrices.") dim = HL[0].shape[-1] # number of basis matrices N = dim**2 - (1 if traceless else 0) if N == 0: return np.empty((0, dim, dim)) # Prepare for differences in sparse and dense algebra if sparse: def basis(): return matrix_basis(dim, traceless=traceless, sparse=True) vstack = scipy.sparse.vstack bmat = scipy.sparse.bmat # Cast it to coo format and reshape def flatten(x): return scipy.sparse.coo_matrix(x).reshape((x.shape[0] * x.shape[1], 1)) else: def basis(): return matrix_basis(dim, traceless=traceless, sparse=False) vstack = np.vstack bmat = np.block def flatten(x): return x.reshape((-1, 1)) # Calculate coefficients from commutators of the basis matrices null_mat = [] for hL, hR, conj in zip(HL, HR, conjugate): if conj: row = [flatten(mat @ hL - hR @ mat.conj()) for mat in basis()] else: row = [flatten(mat @ hL - hR @ mat) for mat in basis()] null_mat.append(row) null_mat = bmat(null_mat) if hermitian: # Simultaneaous equations for real and immaginary parts # Lapack guarantees that the SVD of a real matrix is real null_mat = vstack((null_mat.real, null_mat.imag)) # Find the null space of null_mat ns = nullspace(null_mat, sparse=sparse, k_max=k_max) # Make all solutions # 'ij,jkl->ikl' def basis(): return matrix_basis(dim, traceless=traceless, sparse=False) return np.array([sum((v * mat for v, mat in zip(vec, basis()))) for vec in ns])
[docs] def rref(m, rtol=1e-3, return_S=False): """Bring a matrix to reduced row echelon form. Parameters ---------- m: 2D numpy.ndarray The matrix on which to perform row reduction. rtol: float The tolerance relative to the largest matrix element for what is considered a zero entry. return_S: bool Whether to return the matrix of linear operations that brings the input matrix to reduced row echelon form. Returns --------- red: 2D numpy.ndarray The input matrix in reduced row echelon form. S: 2D numpy.ndarray The matrix of operations that brings the input matrix to reduced row echelon form. Only returned if return_S = True. Note: If r is the reduced row echelon form of the input matrix m, then m = S.dot(r). """ red = np.array(m) rows, columns = red.shape # Define row operations def pivot(i, j): P = np.eye(rows, dtype=complex) P[i, i] = P[j, j] = 0 P[i, j] = P[j, i] = 1 return P def rowsub(i, coeffs): P = np.eye(rows, dtype=complex) P[:, i] -= coeffs P[i, i] = 1 return P def rowmul(i, coeff): P = np.eye(rows, dtype=complex) P[i, i] = coeff return P # Stop if m is identically zero if np.allclose(red, 0): return red S = np.eye(rows, dtype=complex) tol = np.max(np.abs(m)) * rtol lead = 0 for r in range(rows): # find leading value below row r in lead column while lead < columns: i = np.argmax(np.abs(red[r:, lead])) + r lv = red[i, lead] if np.abs(lv) < tol: lead += 1 continue else: break else: break T = pivot(i, r) S = T.dot(S) red = T.dot(red) T = rowmul(r, 1 / lv) S = T.dot(S) red = T.dot(red) T = rowsub(r, red[:, lead]) S = T.dot(S) red = T.dot(red) lead += 1 if return_S: return red, S else: return red
[docs] def sparse_basis(bas, num_digits=3, reals=False): """Reduce the number of nonzero entries in a (basis) matrix by rounding the matrix and bringing it to reduced row echelon form. Parameters ---------- bas: numpy.ndarray The (basis) matrix to prettify and round. num_digits: positive integer Number of significant digits to which to round. reals: boolean If True, the real and imaginary parts of the matrix are treated separately. Returns ------- A rounded, sparsified matrix with the same row span as the input matrix. This is an attempt at matrix sparsification, i.e. we attempt to construct the sparsest matrix that has the same row rank. Note that the reduction to row echelon form is numerically unstable, so this function should be used with caution.""" if len(bas) < 1: return bas if reals: re = np.real(bas) im = np.imag(bas) bas = np.hstack((re, im)) bas_mat = rref(bas, rtol=10 ** (-num_digits)) if reals: re, im = np.split(bas_mat, 2, axis=1) bas_mat = re + 1j * im # Round to num_digits and return only nonvanishing rows bas_mat = np.round(bas_mat, num_digits) bas_mat = np.vstack( [row for row in bas_mat if not np.allclose(row, 0, atol=10 ** (-num_digits))] ) if len(bas_mat) < len(bas): warnings.warn( "Removed linearly dependent terms from the family during sparsification. " + "Resulting family will contain fewer members." ) return np.round(bas_mat, num_digits)
def _inv_int(A): """Invert an integer square matrix A.""" _A = ta.array(A, int) if A == np.empty((0, 0)): return A if _A != A or abs(la.det(A)) != 1: raise ValueError("Input needs to be an invertible integer matrix") return ta.array(np.round(la.inv(_A)), int)