atmos_flux_inversion.linalg module

All the linear algebra details for the package.

I figured it would be more useful to have it separate

class atmos_flux_inversion.linalg.DaskKroneckerProductOperator(*args, **kwargs)[source]

Bases: atmos_flux_inversion.linalg_interface.DaskLinearOperator

Operator for Kronecker product.

Uses the \(O(n^{2.5})\) algorithm of Yadav and Michalak (2013) to make memory and computational requirements practical.

Left argument to Kronecker product must be array.

References

V. Yadav and A.M. Michalak. “Improving computational efficiency in large linear inverse problems: an example from carbon dioxide flux estimation” Geosci. Model Dev. 2013. Vol 6, pp. 583–590. URL: https://www.geosci-model-dev.net/6/583/2013 doi:10.5194/gmd-6-583-2013.

quadratic_form(mat)[source]

Calculate the quadratic form mat.T @ self @ mat.

Parameters

mat (array_like[N, M]) –

Returns

The product mat.T @ self @ mat

Return type

array_like[M, M]

Raises

  • TypeError – if mat is not an array or if self is not square

  • ValueError – if the shapes of self and mat are not compatible

Notes

Implementation depends on Kronecker structure, using the _matmat() algorithm for self @ mat. If mat is a lazy dask array, this implementation will load it multiple times to avoid dask keeping it all in memory.

This function uses dask for the splitting, multiplication, and addition, which defaults to using all available cores.

sqrt()[source]

Find an operator S such that S.T @ S == self.

Requires self be symmetric.

Returns

The square root S of the operator S.T @ S == self

Return type

DaskKroneckerProductOperator

Raises
class atmos_flux_inversion.linalg.DiagonalOperator(*args, **kwargs)[source]

Bases: atmos_flux_inversion.linalg.SelfAdjointLinearOperator

Operator with entries only on the diagonal.

solve(vector)[source]

Solve A @ x == vector.

Parameters

vector (array_like) –

Returns

Solution of self @ x == vec

Return type

array_like

sqrt()[source]

Find S such that S.T @ S == self.

atmos_flux_inversion.linalg.NEAR_ZERO = 1e-20

Where correlations are rounded to zero.

The method of assuring positive definiteness increases some values away from zero due to roundoff. Values that were originally smaller than this are reset to zero.

See also

atmos_flux_inversion.correlations.NEAR_ZERO

atmos_flux_inversion.linalg.OPTIMAL_ELEMENTS = 65536

Maximum size for the approximate square root of a LinearOperator.

Bounding the eigenvector matrix by this should keep everything in memory.

atmos_flux_inversion.linalg.ROUNDOFF = 1e-10

Determines how many terms should be kept from SVD.

class atmos_flux_inversion.linalg.SchmidtKroneckerProduct(*args, **kwargs)[source]

Bases: atmos_flux_inversion.linalg_interface.DaskLinearOperator

Kronecker product of two operators using Schmidt decomposition.

This works best when the input vectors are nearly Kronecker products as well, dominated by some underlying structure with small variations. One example would be average net flux + trend in net flux + average daily cycle + daily cycle timing variations across domain + localized events + …

Multiplications are roughly the same time complexity class as with an explicit Kronecker Product, perhaps a factor of two or three slower in the best case, but the memory requirements are \(N_1^2 + N_2^2\) rather than \((N_1 * N_2)^2\), plus this approach works with sparse matrices and other LinearOperators which can further reduce the memory requirements and may decrease the time complexity.

Forming the Kronecker product from the component vectors currently requires the whole thing to be in memory, so a new implementation of kron would be needed to take advantage of this. There may be some difficulties with the dask cache getting flushed and causing repeat work in this case. I don’t know how to get around this.

class atmos_flux_inversion.linalg.SelfAdjointLinearOperator(*args, **kwargs)[source]

Bases: atmos_flux_inversion.linalg_interface.DaskLinearOperator

Self-adjoint linear operators.

Provides _rmatvec() and _adjoint() methods.

atmos_flux_inversion.linalg.kron(matrix1, matrix2)[source]

Kronecker product of two matrices.

Parameters

  • matrix1 (array_like[M, N]) – One of the matrixes for the product

  • matrix2 (array_like[J, K]) – The other matrix for the product

Returns

The kronecker product of the matrices

Return type

array_like[M*J, N*K]

See also

scipy.linalg.kron()

Where I got the overview of the implementation.

atmos_flux_inversion.linalg.linop_solve(operator, arr)[source]

Solve operator @ x = arr.

deal with arr possibly having multiple columns.

Parameters

  • operator (LinearOperator) –

  • arr (array_like) –

Returns

The solution to the linear equation.

Return type

array_like

atmos_flux_inversion.linalg.matrix_sqrt(mat)[source]

Find a matrix S such that S.T @ S == mat.

Parameters

mat (array_like or LinearOperator) – The square matrix for which the square root is desired

Returns

A matrix such that S.T @ S == mat

Return type

array_like or LinearOperator

Raises
atmos_flux_inversion.linalg.schmidt_decomposition(vector, dim1, dim2)[source]

Decompose a state vector into a sum of Kronecker products.

Parameters

  • vector (array_like[dim1 * dim2]) –

  • dim2 (dim1,) –

Returns

The rows form the separate vectors. The weights are guaranteed to be greater than zero

Return type

tuple of (weights, unit_vecs[dim1], unit_vecs[dim2]

Raises

ValueError – if vector isn’t actually a vector.

Notes

Algorithm from stackexchange: https://physics.stackexchange.com/questions/251522/how-do-you-find-a-schmidt-basis-and-how-can-the-schmidt-decomposition-be-used-f Also from Mathematica code I wrote based on description in the green Quantum Computation book in the reading library

atmos_flux_inversion.linalg.solve(arr1, arr2)[source]

Solve arr1 @ x = arr2 for x.

Parameters

  • arr1 (array_like[N, N]) – A square matrix

  • arr2 (array_like[N]) – A vector

Returns

The solution to the linear equation

Return type

array_like[N]

Raises

  • ValueError – if the dimensions do not match up

  • LinAlgError – if arr1 is not square