pyerrors.linalg
1import numpy as np 2import autograd.numpy as anp # Thinly-wrapped numpy 3from .obs import derived_observable, CObs, Obs, import_jackknife 4 5 6def matmul(*operands): 7 """Matrix multiply all operands. 8 9 Parameters 10 ---------- 11 operands : numpy.ndarray 12 Arbitrary number of 2d-numpy arrays which can be real or complex 13 Obs valued. 14 15 This implementation is faster compared to standard multiplication via the @ operator. 16 """ 17 if any(isinstance(o[0, 0], CObs) for o in operands): 18 extended_operands = [] 19 for op in operands: 20 tmp = np.vectorize(lambda x: (np.real(x), np.imag(x)))(op) 21 extended_operands.append(tmp[0]) 22 extended_operands.append(tmp[1]) 23 24 def multi_dot(operands, part): 25 stack_r = operands[0] 26 stack_i = operands[1] 27 for op_r, op_i in zip(operands[2::2], operands[3::2]): 28 tmp_r = stack_r @ op_r - stack_i @ op_i 29 tmp_i = stack_r @ op_i + stack_i @ op_r 30 31 stack_r = tmp_r 32 stack_i = tmp_i 33 34 if part == 'Real': 35 return stack_r 36 else: 37 return stack_i 38 39 def multi_dot_r(operands): 40 return multi_dot(operands, 'Real') 41 42 def multi_dot_i(operands): 43 return multi_dot(operands, 'Imag') 44 45 Nr = derived_observable(multi_dot_r, extended_operands, array_mode=True) 46 Ni = derived_observable(multi_dot_i, extended_operands, array_mode=True) 47 48 res = np.empty_like(Nr) 49 for (n, m), entry in np.ndenumerate(Nr): 50 res[n, m] = CObs(Nr[n, m], Ni[n, m]) 51 52 return res 53 else: 54 def multi_dot(operands): 55 stack = operands[0] 56 for op in operands[1:]: 57 stack = stack @ op 58 return stack 59 return derived_observable(multi_dot, operands, array_mode=True) 60 61 62def jack_matmul(*operands): 63 """Matrix multiply both operands making use of the jackknife approximation. 64 65 Parameters 66 ---------- 67 operands : numpy.ndarray 68 Arbitrary number of 2d-numpy arrays which can be real or complex 69 Obs valued. 70 71 For large matrices this is considerably faster compared to matmul. 72 """ 73 74 def _exp_to_jack(matrix): 75 base_matrix = np.empty_like(matrix) 76 for index, entry in np.ndenumerate(matrix): 77 base_matrix[index] = entry.export_jackknife() 78 return base_matrix 79 80 def _imp_from_jack(matrix, name, idl): 81 base_matrix = np.empty_like(matrix) 82 for index, entry in np.ndenumerate(matrix): 83 base_matrix[index] = import_jackknife(entry, name, [idl]) 84 return base_matrix 85 86 def _exp_to_jack_c(matrix): 87 base_matrix = np.empty_like(matrix) 88 for index, entry in np.ndenumerate(matrix): 89 base_matrix[index] = entry.real.export_jackknife() + 1j * entry.imag.export_jackknife() 90 return base_matrix 91 92 def _imp_from_jack_c(matrix, name, idl): 93 base_matrix = np.empty_like(matrix) 94 for index, entry in np.ndenumerate(matrix): 95 base_matrix[index] = CObs(import_jackknife(entry.real, name, [idl]), 96 import_jackknife(entry.imag, name, [idl])) 97 return base_matrix 98 99 if any(isinstance(o.flat[0], CObs) for o in operands): 100 name = operands[0].flat[0].real.names[0] 101 idl = operands[0].flat[0].real.idl[name] 102 103 r = _exp_to_jack_c(operands[0]) 104 for op in operands[1:]: 105 if isinstance(op.flat[0], CObs): 106 r = r @ _exp_to_jack_c(op) 107 else: 108 r = r @ op 109 return _imp_from_jack_c(r, name, idl) 110 else: 111 name = operands[0].flat[0].names[0] 112 idl = operands[0].flat[0].idl[name] 113 114 r = _exp_to_jack(operands[0]) 115 for op in operands[1:]: 116 if isinstance(op.flat[0], Obs): 117 r = r @ _exp_to_jack(op) 118 else: 119 r = r @ op 120 return _imp_from_jack(r, name, idl) 121 122 123def einsum(subscripts, *operands): 124 """Wrapper for numpy.einsum 125 126 Parameters 127 ---------- 128 subscripts : str 129 Subscripts for summation (see numpy documentation for details) 130 operands : numpy.ndarray 131 Arbitrary number of 2d-numpy arrays which can be real or complex 132 Obs valued. 133 """ 134 135 def _exp_to_jack(matrix): 136 base_matrix = [] 137 for index, entry in np.ndenumerate(matrix): 138 base_matrix.append(entry.export_jackknife()) 139 return np.asarray(base_matrix).reshape(matrix.shape + base_matrix[0].shape) 140 141 def _exp_to_jack_c(matrix): 142 base_matrix = [] 143 for index, entry in np.ndenumerate(matrix): 144 base_matrix.append(entry.real.export_jackknife() + 1j * entry.imag.export_jackknife()) 145 return np.asarray(base_matrix).reshape(matrix.shape + base_matrix[0].shape) 146 147 def _imp_from_jack(matrix, name, idl): 148 base_matrix = np.empty(shape=matrix.shape[:-1], dtype=object) 149 for index in np.ndindex(matrix.shape[:-1]): 150 base_matrix[index] = import_jackknife(matrix[index], name, [idl]) 151 return base_matrix 152 153 def _imp_from_jack_c(matrix, name, idl): 154 base_matrix = np.empty(shape=matrix.shape[:-1], dtype=object) 155 for index in np.ndindex(matrix.shape[:-1]): 156 base_matrix[index] = CObs(import_jackknife(matrix[index].real, name, [idl]), 157 import_jackknife(matrix[index].imag, name, [idl])) 158 return base_matrix 159 160 for op in operands: 161 if isinstance(op.flat[0], CObs): 162 name = op.flat[0].real.names[0] 163 idl = op.flat[0].real.idl[name] 164 break 165 elif isinstance(op.flat[0], Obs): 166 name = op.flat[0].names[0] 167 idl = op.flat[0].idl[name] 168 break 169 170 conv_operands = [] 171 for op in operands: 172 if isinstance(op.flat[0], CObs): 173 conv_operands.append(_exp_to_jack_c(op)) 174 elif isinstance(op.flat[0], Obs): 175 conv_operands.append(_exp_to_jack(op)) 176 else: 177 conv_operands.append(op) 178 179 tmp_subscripts = ','.join([o + '...' for o in subscripts.split(',')]) 180 extended_subscripts = '->'.join([o + '...' for o in tmp_subscripts.split('->')[:-1]] + [tmp_subscripts.split('->')[-1]]) 181 einsum_path = np.einsum_path(extended_subscripts, *conv_operands, optimize='optimal')[0] 182 jack_einsum = np.einsum(extended_subscripts, *conv_operands, optimize=einsum_path) 183 184 if jack_einsum.dtype == complex: 185 result = _imp_from_jack_c(jack_einsum, name, idl) 186 elif jack_einsum.dtype == float: 187 result = _imp_from_jack(jack_einsum, name, idl) 188 else: 189 raise Exception("Result has unexpected datatype") 190 191 if result.shape == (): 192 return result.flat[0] 193 else: 194 return result 195 196 197def inv(x): 198 """Inverse of Obs or CObs valued matrices.""" 199 return _mat_mat_op(anp.linalg.inv, x) 200 201 202def cholesky(x): 203 """Cholesky decomposition of Obs valued matrices.""" 204 if any(isinstance(o, CObs) for o in x.ravel()): 205 raise Exception("Cholesky decomposition is not implemented for CObs.") 206 return _mat_mat_op(anp.linalg.cholesky, x) 207 208 209def det(x): 210 """Determinant of Obs valued matrices.""" 211 return _scalar_mat_op(anp.linalg.det, x) 212 213 214def _scalar_mat_op(op, obs, **kwargs): 215 """Computes the matrix to scalar operation op to a given matrix of Obs.""" 216 def _mat(x, **kwargs): 217 dim = int(np.sqrt(len(x))) 218 219 mat = [] 220 for i in range(dim): 221 row = [] 222 for j in range(dim): 223 row.append(x[j + dim * i]) 224 mat.append(row) 225 226 return op(anp.array(mat)) 227 228 if isinstance(obs, np.ndarray): 229 raveled_obs = (1 * (obs.ravel())).tolist() 230 else: 231 raise TypeError('Unproper type of input.') 232 return derived_observable(_mat, raveled_obs, **kwargs) 233 234 235def _mat_mat_op(op, obs, **kwargs): 236 """Computes the matrix to matrix operation op to a given matrix of Obs.""" 237 # Use real representation to calculate matrix operations for complex matrices 238 if any(isinstance(o, CObs) for o in obs.ravel()): 239 A = np.empty_like(obs) 240 B = np.empty_like(obs) 241 for (n, m), entry in np.ndenumerate(obs): 242 if hasattr(entry, 'real') and hasattr(entry, 'imag'): 243 A[n, m] = entry.real 244 B[n, m] = entry.imag 245 else: 246 A[n, m] = entry 247 B[n, m] = 0.0 248 big_matrix = np.block([[A, -B], [B, A]]) 249 op_big_matrix = derived_observable(lambda x, **kwargs: op(x), [big_matrix], array_mode=True)[0] 250 dim = op_big_matrix.shape[0] 251 op_A = op_big_matrix[0: dim // 2, 0: dim // 2] 252 op_B = op_big_matrix[dim // 2:, 0: dim // 2] 253 res = np.empty_like(op_A) 254 for (n, m), entry in np.ndenumerate(op_A): 255 res[n, m] = CObs(op_A[n, m], op_B[n, m]) 256 return res 257 else: 258 return derived_observable(lambda x, **kwargs: op(x), [obs], array_mode=True)[0] 259 260 261def eigh(obs, **kwargs): 262 """Computes the eigenvalues and eigenvectors of a given hermitian matrix of Obs according to np.linalg.eigh.""" 263 w = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[0], obs) 264 v = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[1], obs) 265 return w, v 266 267 268def eig(obs, **kwargs): 269 """Computes the eigenvalues of a given matrix of Obs according to np.linalg.eig.""" 270 w = derived_observable(lambda x, **kwargs: anp.real(anp.linalg.eig(x)[0]), obs) 271 return w 272 273 274def pinv(obs, **kwargs): 275 """Computes the Moore-Penrose pseudoinverse of a matrix of Obs.""" 276 return derived_observable(lambda x, **kwargs: anp.linalg.pinv(x), obs) 277 278 279def svd(obs, **kwargs): 280 """Computes the singular value decomposition of a matrix of Obs.""" 281 u = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[0], obs) 282 s = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[1], obs) 283 vh = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[2], obs) 284 return (u, s, vh)
def
matmul(*operands):
7def matmul(*operands): 8 """Matrix multiply all operands. 9 10 Parameters 11 ---------- 12 operands : numpy.ndarray 13 Arbitrary number of 2d-numpy arrays which can be real or complex 14 Obs valued. 15 16 This implementation is faster compared to standard multiplication via the @ operator. 17 """ 18 if any(isinstance(o[0, 0], CObs) for o in operands): 19 extended_operands = [] 20 for op in operands: 21 tmp = np.vectorize(lambda x: (np.real(x), np.imag(x)))(op) 22 extended_operands.append(tmp[0]) 23 extended_operands.append(tmp[1]) 24 25 def multi_dot(operands, part): 26 stack_r = operands[0] 27 stack_i = operands[1] 28 for op_r, op_i in zip(operands[2::2], operands[3::2]): 29 tmp_r = stack_r @ op_r - stack_i @ op_i 30 tmp_i = stack_r @ op_i + stack_i @ op_r 31 32 stack_r = tmp_r 33 stack_i = tmp_i 34 35 if part == 'Real': 36 return stack_r 37 else: 38 return stack_i 39 40 def multi_dot_r(operands): 41 return multi_dot(operands, 'Real') 42 43 def multi_dot_i(operands): 44 return multi_dot(operands, 'Imag') 45 46 Nr = derived_observable(multi_dot_r, extended_operands, array_mode=True) 47 Ni = derived_observable(multi_dot_i, extended_operands, array_mode=True) 48 49 res = np.empty_like(Nr) 50 for (n, m), entry in np.ndenumerate(Nr): 51 res[n, m] = CObs(Nr[n, m], Ni[n, m]) 52 53 return res 54 else: 55 def multi_dot(operands): 56 stack = operands[0] 57 for op in operands[1:]: 58 stack = stack @ op 59 return stack 60 return derived_observable(multi_dot, operands, array_mode=True)
Matrix multiply all operands.
Parameters
- operands (numpy.ndarray): Arbitrary number of 2d-numpy arrays which can be real or complex Obs valued.
- This implementation is faster compared to standard multiplication via the @ operator.
def
jack_matmul(*operands):
63def jack_matmul(*operands): 64 """Matrix multiply both operands making use of the jackknife approximation. 65 66 Parameters 67 ---------- 68 operands : numpy.ndarray 69 Arbitrary number of 2d-numpy arrays which can be real or complex 70 Obs valued. 71 72 For large matrices this is considerably faster compared to matmul. 73 """ 74 75 def _exp_to_jack(matrix): 76 base_matrix = np.empty_like(matrix) 77 for index, entry in np.ndenumerate(matrix): 78 base_matrix[index] = entry.export_jackknife() 79 return base_matrix 80 81 def _imp_from_jack(matrix, name, idl): 82 base_matrix = np.empty_like(matrix) 83 for index, entry in np.ndenumerate(matrix): 84 base_matrix[index] = import_jackknife(entry, name, [idl]) 85 return base_matrix 86 87 def _exp_to_jack_c(matrix): 88 base_matrix = np.empty_like(matrix) 89 for index, entry in np.ndenumerate(matrix): 90 base_matrix[index] = entry.real.export_jackknife() + 1j * entry.imag.export_jackknife() 91 return base_matrix 92 93 def _imp_from_jack_c(matrix, name, idl): 94 base_matrix = np.empty_like(matrix) 95 for index, entry in np.ndenumerate(matrix): 96 base_matrix[index] = CObs(import_jackknife(entry.real, name, [idl]), 97 import_jackknife(entry.imag, name, [idl])) 98 return base_matrix 99 100 if any(isinstance(o.flat[0], CObs) for o in operands): 101 name = operands[0].flat[0].real.names[0] 102 idl = operands[0].flat[0].real.idl[name] 103 104 r = _exp_to_jack_c(operands[0]) 105 for op in operands[1:]: 106 if isinstance(op.flat[0], CObs): 107 r = r @ _exp_to_jack_c(op) 108 else: 109 r = r @ op 110 return _imp_from_jack_c(r, name, idl) 111 else: 112 name = operands[0].flat[0].names[0] 113 idl = operands[0].flat[0].idl[name] 114 115 r = _exp_to_jack(operands[0]) 116 for op in operands[1:]: 117 if isinstance(op.flat[0], Obs): 118 r = r @ _exp_to_jack(op) 119 else: 120 r = r @ op 121 return _imp_from_jack(r, name, idl)
Matrix multiply both operands making use of the jackknife approximation.
Parameters
- operands (numpy.ndarray): Arbitrary number of 2d-numpy arrays which can be real or complex Obs valued.
- For large matrices this is considerably faster compared to matmul.
def
einsum(subscripts, *operands):
124def einsum(subscripts, *operands): 125 """Wrapper for numpy.einsum 126 127 Parameters 128 ---------- 129 subscripts : str 130 Subscripts for summation (see numpy documentation for details) 131 operands : numpy.ndarray 132 Arbitrary number of 2d-numpy arrays which can be real or complex 133 Obs valued. 134 """ 135 136 def _exp_to_jack(matrix): 137 base_matrix = [] 138 for index, entry in np.ndenumerate(matrix): 139 base_matrix.append(entry.export_jackknife()) 140 return np.asarray(base_matrix).reshape(matrix.shape + base_matrix[0].shape) 141 142 def _exp_to_jack_c(matrix): 143 base_matrix = [] 144 for index, entry in np.ndenumerate(matrix): 145 base_matrix.append(entry.real.export_jackknife() + 1j * entry.imag.export_jackknife()) 146 return np.asarray(base_matrix).reshape(matrix.shape + base_matrix[0].shape) 147 148 def _imp_from_jack(matrix, name, idl): 149 base_matrix = np.empty(shape=matrix.shape[:-1], dtype=object) 150 for index in np.ndindex(matrix.shape[:-1]): 151 base_matrix[index] = import_jackknife(matrix[index], name, [idl]) 152 return base_matrix 153 154 def _imp_from_jack_c(matrix, name, idl): 155 base_matrix = np.empty(shape=matrix.shape[:-1], dtype=object) 156 for index in np.ndindex(matrix.shape[:-1]): 157 base_matrix[index] = CObs(import_jackknife(matrix[index].real, name, [idl]), 158 import_jackknife(matrix[index].imag, name, [idl])) 159 return base_matrix 160 161 for op in operands: 162 if isinstance(op.flat[0], CObs): 163 name = op.flat[0].real.names[0] 164 idl = op.flat[0].real.idl[name] 165 break 166 elif isinstance(op.flat[0], Obs): 167 name = op.flat[0].names[0] 168 idl = op.flat[0].idl[name] 169 break 170 171 conv_operands = [] 172 for op in operands: 173 if isinstance(op.flat[0], CObs): 174 conv_operands.append(_exp_to_jack_c(op)) 175 elif isinstance(op.flat[0], Obs): 176 conv_operands.append(_exp_to_jack(op)) 177 else: 178 conv_operands.append(op) 179 180 tmp_subscripts = ','.join([o + '...' for o in subscripts.split(',')]) 181 extended_subscripts = '->'.join([o + '...' for o in tmp_subscripts.split('->')[:-1]] + [tmp_subscripts.split('->')[-1]]) 182 einsum_path = np.einsum_path(extended_subscripts, *conv_operands, optimize='optimal')[0] 183 jack_einsum = np.einsum(extended_subscripts, *conv_operands, optimize=einsum_path) 184 185 if jack_einsum.dtype == complex: 186 result = _imp_from_jack_c(jack_einsum, name, idl) 187 elif jack_einsum.dtype == float: 188 result = _imp_from_jack(jack_einsum, name, idl) 189 else: 190 raise Exception("Result has unexpected datatype") 191 192 if result.shape == (): 193 return result.flat[0] 194 else: 195 return result
Wrapper for numpy.einsum
Parameters
- subscripts (str): Subscripts for summation (see numpy documentation for details)
- operands (numpy.ndarray): Arbitrary number of 2d-numpy arrays which can be real or complex Obs valued.
def
inv(x):
198def inv(x): 199 """Inverse of Obs or CObs valued matrices.""" 200 return _mat_mat_op(anp.linalg.inv, x)
Inverse of Obs or CObs valued matrices.
def
cholesky(x):
203def cholesky(x): 204 """Cholesky decomposition of Obs valued matrices.""" 205 if any(isinstance(o, CObs) for o in x.ravel()): 206 raise Exception("Cholesky decomposition is not implemented for CObs.") 207 return _mat_mat_op(anp.linalg.cholesky, x)
Cholesky decomposition of Obs valued matrices.
def
det(x):
210def det(x): 211 """Determinant of Obs valued matrices.""" 212 return _scalar_mat_op(anp.linalg.det, x)
Determinant of Obs valued matrices.
def
eigh(obs, **kwargs):
262def eigh(obs, **kwargs): 263 """Computes the eigenvalues and eigenvectors of a given hermitian matrix of Obs according to np.linalg.eigh.""" 264 w = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[0], obs) 265 v = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[1], obs) 266 return w, v
Computes the eigenvalues and eigenvectors of a given hermitian matrix of Obs according to np.linalg.eigh.
def
eig(obs, **kwargs):
269def eig(obs, **kwargs): 270 """Computes the eigenvalues of a given matrix of Obs according to np.linalg.eig.""" 271 w = derived_observable(lambda x, **kwargs: anp.real(anp.linalg.eig(x)[0]), obs) 272 return w
Computes the eigenvalues of a given matrix of Obs according to np.linalg.eig.
def
pinv(obs, **kwargs):
275def pinv(obs, **kwargs): 276 """Computes the Moore-Penrose pseudoinverse of a matrix of Obs.""" 277 return derived_observable(lambda x, **kwargs: anp.linalg.pinv(x), obs)
Computes the Moore-Penrose pseudoinverse of a matrix of Obs.
def
svd(obs, **kwargs):
280def svd(obs, **kwargs): 281 """Computes the singular value decomposition of a matrix of Obs.""" 282 u = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[0], obs) 283 s = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[1], obs) 284 vh = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[2], obs) 285 return (u, s, vh)
Computes the singular value decomposition of a matrix of Obs.