pyerrors.linalg

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

 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)

 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)

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

197def inv(x):
198    """Inverse of Obs or CObs valued matrices."""
199    return _mat_mat_op(anp.linalg.inv, x)

202def cholesky(x):
203    """Cholesky decomposition of Obs or CObs valued matrices."""
204    return _mat_mat_op(anp.linalg.cholesky, x)

207def det(x):
208    """Determinant of Obs valued matrices."""
209    return _scalar_mat_op(anp.linalg.det, x)

259def eigh(obs, **kwargs):
260    """Computes the eigenvalues and eigenvectors of a given hermitian matrix of Obs according to np.linalg.eigh."""
261    w = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[0], obs)
262    v = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[1], obs)
263    return w, v

266def eig(obs, **kwargs):
267    """Computes the eigenvalues of a given matrix of Obs according to np.linalg.eig."""
268    w = derived_observable(lambda x, **kwargs: anp.real(anp.linalg.eig(x)[0]), obs)
269    return w

272def pinv(obs, **kwargs):
273    """Computes the Moore-Penrose pseudoinverse of a matrix of Obs."""
274    return derived_observable(lambda x, **kwargs: anp.linalg.pinv(x), obs)

277def svd(obs, **kwargs):
278    """Computes the singular value decomposition of a matrix of Obs."""
279    u = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[0], obs)
280    s = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[1], obs)
281    vh = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[2], obs)
282    return (u, s, vh)