solvers package#

Submodules#

PyHydroGeophysX.solvers.linear_solvers module#

Linear solvers for geophysical inversion.

class PyHydroGeophysX.solvers.linear_solvers.CGLSSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]#

Bases: LinearSolver

CGLS (Conjugate Gradient Least Squares) solver.

class PyHydroGeophysX.solvers.linear_solvers.IterativeRefinement(max_iterations=5, tolerance=1e-10, use_double_precision=True)[source]#

Bases: object

Iterative refinement to improve accuracy of a solution to a linear system.

refine(A, b, x0, solver_func)[source]#

Perform iterative refinement.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • x0 – Initial solution

  • solver_func – Function that solves A*x = b

Returns:

Improved solution

class PyHydroGeophysX.solvers.linear_solvers.LSQRSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]#

Bases: LinearSolver

LSQR solver for least squares problems.

class PyHydroGeophysX.solvers.linear_solvers.LinearSolver(method='cgls', max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]#

Bases: object

Base class for linear system solvers.

solve(A, b, x0=None)[source]#

Solve linear system Ax = b.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • x0 – Initial guess (None for zeros)

Returns:

Solution vector

class PyHydroGeophysX.solvers.linear_solvers.RRLSQRSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.1, verbose=False)[source]#

Bases: LinearSolver

Regularized LSQR solver.

class PyHydroGeophysX.solvers.linear_solvers.RRLSSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]#

Bases: LinearSolver

Range-Restricted Least Squares solver.

class PyHydroGeophysX.solvers.linear_solvers.StackedSystem(*args: Any, **kwargs: Any)[source]#

Bases: LinearOperator

A stacked least-squares system [B_1; B_2; ...] held as its row blocks.

An inversion’s Gauss-Newton system, [W_d J; sqrt(lambda) W_m; ...], is a dense Jacobian over a sparse regularization. Stacked into one dense array the regularization rows cost C x n doubles - 1.5 n^2 for a pyGIMLi smoothness operator, 2.7 GB at n = 15000 - and every product with the system multiplied their zeros. Held as blocks, each block keeps the form it was built in.

The least-squares solvers only multiply by the system and its transpose, which this does as a LinearOperator. normal_matrix gives an SPD solver A^T A summed block by block, columns keeps a subset of the columns (cells held fixed), and toarray assembles the dense array for a solver that needs one, such as the GPU path.

column_norms() → ndarray[source]#

The Euclidean norm of every column.

columns(keep: Any) → StackedSystem[source]#

The system restricted to the columns keep selects.

normal_matrix() → ndarray[source]#

A^T A as a dense array, summed block by block.

toarray() → ndarray[source]#
class PyHydroGeophysX.solvers.linear_solvers.TikhonvRegularization(regularization_matrix=None, alpha=1.0, regularization_type='identity')[source]#

Bases: object

Tikhonov regularization for ill-posed problems.

apply(A, b, solver=None)[source]#

Apply Tikhonov regularization to the linear system.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • solver – Solver to use (None for direct solver)

Returns:

Regularized solution

create_regularization_matrix(n)[source]#

Create regularization matrix based on the selected type.

Parameters:

n – Size of model vector

Returns:

Regularization matrix

PyHydroGeophysX.solvers.linear_solvers.block_tridiagonal_cholesky_solve(diagonal: Any, lower: Any, b: Any, *, overwrite: bool = False) → Any[source]#

Solve H x = b for a symmetric positive definite block-tridiagonal H.

diagonal[k] is the dense block H[k, k]; lower[k] is the block H[k + 1, k] below it, either a dense matrix or a vector standing for that diagonal matrix - the form in which a temporal-difference penalty couples adjacent surveys of a time-lapse inversion.

The Cholesky factor of such a matrix is block lower-bidiagonal, so this is the dense Cholesky factorization with its zero blocks left out: a Cholesky of each Schur complement, a triangular solve for the block below it and a symmetric rank update of the next diagonal block - the operations a blocked dense factorization performs on the blocks that are not zero, in O(N n^3) work and O(N n^2) memory rather than O((N n)^3) and O((N n)^2).

With overwrite the factors are written into the diagonal blocks (float64, contiguous) and nothing of that size is allocated. A block that is not positive definite raises LinAlgError: as for spd_cholesky with overwrite_a=True, there is no ridge retry.

>>> H = np.array([[4.0, 1.0, -0.5, 0.0], [1.0, 3.0, 0.0, -0.5],
...               [-0.5, 0.0, 5.0, 2.0], [0.0, -0.5, 2.0, 6.0]])
>>> x = block_tridiagonal_cholesky_solve(
...     [H[:2, :2], H[2:, 2:]], [np.array([-0.5, -0.5])], np.ones(4))
>>> bool(np.allclose(H @ x.ravel(), 1.0))
True
PyHydroGeophysX.solvers.linear_solvers.direct_solver(A: Any, b: Any, method: Any = 'lu', **kwargs: Any) → Any[source]#

Solve a linear system using direct methods.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • method – Direct solver method (‘lu’, ‘qr’, ‘svd’, ‘cholesky’)

  • **kwargs – Additional parameters for specific methods

Returns:

Solution vector

PyHydroGeophysX.solvers.linear_solvers.generalized_solver(A: Any, b: Any, method: Any = 'cgls', x: Any = None, maxiter: Any = 200, tol: Any = 1e-08, verbose: Any = False, damp: Any = 0.0, use_gpu: Any = False, parallel: Any = False, n_jobs: Any = -1, *, overwrite_a: Any = False) → Any[source]#

Generalized solver for Ax = b with optional GPU acceleration and parallelism.

Parameters:#

Aarray_like or sparse matrix

The system matrix. Two shapes of problem reach this function, and they need opposite solver families:

  1. A stacked least-squares system, [W_d J; sqrt(lambda) W_m; ...], with a matching stacked right-hand side. Use a least-squares method.

  2. A Gauss-Newton normal matrix, J^T W_d^T W_d J + lambda W_m^T W_m + ..., which is square and symmetric (positive definite when the combined operators constrain every model direction), paired with a gradient right-hand side -g. Use an SPD method. A least-squares method applied to this case solves A^T A x = A^T b instead, which squares the condition number and spends two matrix-vector products per iteration where one would do; passing one is warned about once per process.

barray_like

Right-hand side vector.

methodstr, optional

Solver method. Default is ‘cgls’. Iterative: ‘lsqr’, ‘rrlsqr’, ‘cgls’, ‘rrls’ SciPy: ‘scipy_lsqr’, ‘scipy_lsmr’, ‘precond_lsmr’ SPD, CPU-only, for a square symmetric positive definite A: ‘spd_cholesky’ (direct factorization), ‘spd_cg’ (iterative, for a matrix too large to factor)

xarray_like, optional

Initial guess for the solution. If None, zeros are used.

maxiterint, optional

Maximum number of iterations.

tolfloat, optional

Convergence tolerance.

verbosebool, optional

Print progress information every 10 iterations.

dampfloat, optional

Damping factor (Tikhonov regularization).

use_gpubool, optional

Use GPU acceleration with CuPy (if available).

parallelbool, optional

Use parallel CPU computations.

n_jobsint, optional

Number of parallel jobs (if parallel is True).

overwrite_abool, optional

Keyword-only. ‘spd_cholesky’ only. Factor in the caller’s own buffer, allowing compatible input storage to be overwritten. SciPy may still copy or allocate workspace. Set it only when A is not needed after the call. Ignored by every other method.

Returns:#

xarray_like

The computed solution vector.

PyHydroGeophysX.solvers.linear_solvers.get_optimal_solver(A: Any, b: Any, estimate_condition: Any = True, time_limit: Any = None, memory_limit: Any = None) → Any[source]#

Automatically select the optimal solver for a given linear system.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • estimate_condition – Whether to estimate condition number

  • time_limit – Maximum allowed solution time (seconds)

  • memory_limit – Maximum allowed memory usage (bytes)

Returns:

Tuple of (solver_object, solver_info)

PyHydroGeophysX.solvers.linear_solvers.spd_solve(A: Any, B: Any, what: Any = 'matrix') → Any[source]#

Solve A X = B for a symmetric positive definite A.

Use this wherever the code would otherwise write inv(A) @ B or pinv(A) @ B. Three reasons it is the better default for a covariance or a Gauss-Newton normal matrix:

  • Cholesky solves the positive-definite system without constructing an explicit SVD pseudoinverse. Runtime depends on size and BLAS backend.

  • pinv discards singular values below its rcond without saying so. On a normal matrix, whose condition number is already the square of the Jacobian’s, that cut can land inside the real spectrum.

  • Finite-precision inverse construction can introduce asymmetry or small negative eigenvalues in ill-conditioned covariance calculations.

A matrix that is genuinely not positive definite falls back to pinv with a warning naming what, so a caller learns which quantity degraded rather than silently receiving a worse answer.

PyHydroGeophysX.solvers.linear_solvers.symmetrize(A: Any) → Any[source]#

Return (A + A.T) / 2.

A covariance or resolution matrix is symmetric in exact arithmetic. Products of the form J^T C J accumulate asymmetry at the rounding level, and np.random.Generator.multivariate_normal rejects the result outright, so the halves are averaged before the matrix leaves the function that built it.

PyHydroGeophysX.solvers.solver module#

Deprecated in 0.5.0 and removed in 0.6.0: a compatibility path whose generalized_solver warns and calls PyHydroGeophysX.solvers.linear_solvers.generalized_solver(), documented above.

Module contents#

Solver exports for inversion and regularized linear-system utilities.

class PyHydroGeophysX.solvers.CGLSSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]#

Bases: LinearSolver

CGLS (Conjugate Gradient Least Squares) solver.

class PyHydroGeophysX.solvers.IterativeRefinement(max_iterations=5, tolerance=1e-10, use_double_precision=True)[source]#

Bases: object

Iterative refinement to improve accuracy of a solution to a linear system.

refine(A, b, x0, solver_func)[source]#

Perform iterative refinement.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • x0 – Initial solution

  • solver_func – Function that solves A*x = b

Returns:

Improved solution

class PyHydroGeophysX.solvers.LSQRSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]#

Bases: LinearSolver

LSQR solver for least squares problems.

class PyHydroGeophysX.solvers.LinearSolver(method='cgls', max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]#

Bases: object

Base class for linear system solvers.

solve(A, b, x0=None)[source]#

Solve linear system Ax = b.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • x0 – Initial guess (None for zeros)

Returns:

Solution vector

class PyHydroGeophysX.solvers.RRLSQRSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.1, verbose=False)[source]#

Bases: LinearSolver

Regularized LSQR solver.

class PyHydroGeophysX.solvers.RRLSSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]#

Bases: LinearSolver

Range-Restricted Least Squares solver.

class PyHydroGeophysX.solvers.TikhonvRegularization(regularization_matrix=None, alpha=1.0, regularization_type='identity')[source]#

Bases: object

Tikhonov regularization for ill-posed problems.

apply(A, b, solver=None)[source]#

Apply Tikhonov regularization to the linear system.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • solver – Solver to use (None for direct solver)

Returns:

Regularized solution

create_regularization_matrix(n)[source]#

Create regularization matrix based on the selected type.

Parameters:

n – Size of model vector

Returns:

Regularization matrix

PyHydroGeophysX.solvers.direct_solver(A: Any, b: Any, method: Any = 'lu', **kwargs: Any) → Any[source]#

Solve a linear system using direct methods.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • method – Direct solver method (‘lu’, ‘qr’, ‘svd’, ‘cholesky’)

  • **kwargs – Additional parameters for specific methods

Returns:

Solution vector

PyHydroGeophysX.solvers.generalized_solver(A: Any, b: Any, method: Any = 'cgls', x: Any = None, maxiter: Any = 200, tol: Any = 1e-08, verbose: Any = False, damp: Any = 0.0, use_gpu: Any = False, parallel: Any = False, n_jobs: Any = -1, *, overwrite_a: Any = False) → Any[source]#

Generalized solver for Ax = b with optional GPU acceleration and parallelism.

Parameters:#

Aarray_like or sparse matrix

The system matrix. Two shapes of problem reach this function, and they need opposite solver families:

  1. A stacked least-squares system, [W_d J; sqrt(lambda) W_m; ...], with a matching stacked right-hand side. Use a least-squares method.

  2. A Gauss-Newton normal matrix, J^T W_d^T W_d J + lambda W_m^T W_m + ..., which is square and symmetric (positive definite when the combined operators constrain every model direction), paired with a gradient right-hand side -g. Use an SPD method. A least-squares method applied to this case solves A^T A x = A^T b instead, which squares the condition number and spends two matrix-vector products per iteration where one would do; passing one is warned about once per process.

barray_like

Right-hand side vector.

methodstr, optional

Solver method. Default is ‘cgls’. Iterative: ‘lsqr’, ‘rrlsqr’, ‘cgls’, ‘rrls’ SciPy: ‘scipy_lsqr’, ‘scipy_lsmr’, ‘precond_lsmr’ SPD, CPU-only, for a square symmetric positive definite A: ‘spd_cholesky’ (direct factorization), ‘spd_cg’ (iterative, for a matrix too large to factor)

xarray_like, optional

Initial guess for the solution. If None, zeros are used.

maxiterint, optional

Maximum number of iterations.

tolfloat, optional

Convergence tolerance.

verbosebool, optional

Print progress information every 10 iterations.

dampfloat, optional

Damping factor (Tikhonov regularization).

use_gpubool, optional

Use GPU acceleration with CuPy (if available).

parallelbool, optional

Use parallel CPU computations.

n_jobsint, optional

Number of parallel jobs (if parallel is True).

overwrite_abool, optional

Keyword-only. ‘spd_cholesky’ only. Factor in the caller’s own buffer, allowing compatible input storage to be overwritten. SciPy may still copy or allocate workspace. Set it only when A is not needed after the call. Ignored by every other method.

Returns:#

xarray_like

The computed solution vector.

PyHydroGeophysX.solvers.get_optimal_solver(A: Any, b: Any, estimate_condition: Any = True, time_limit: Any = None, memory_limit: Any = None) → Any[source]#

Automatically select the optimal solver for a given linear system.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • estimate_condition – Whether to estimate condition number

  • time_limit – Maximum allowed solution time (seconds)

  • memory_limit – Maximum allowed memory usage (bytes)

Returns:

Tuple of (solver_object, solver_info)

PyHydroGeophysX.solvers.spd_solve(A: Any, B: Any, what: Any = 'matrix') → Any[source]#

Solve A X = B for a symmetric positive definite A.

Use this wherever the code would otherwise write inv(A) @ B or pinv(A) @ B. Three reasons it is the better default for a covariance or a Gauss-Newton normal matrix:

  • Cholesky solves the positive-definite system without constructing an explicit SVD pseudoinverse. Runtime depends on size and BLAS backend.

  • pinv discards singular values below its rcond without saying so. On a normal matrix, whose condition number is already the square of the Jacobian’s, that cut can land inside the real spectrum.

  • Finite-precision inverse construction can introduce asymmetry or small negative eigenvalues in ill-conditioned covariance calculations.

A matrix that is genuinely not positive definite falls back to pinv with a warning naming what, so a caller learns which quantity degraded rather than silently receiving a worse answer.

PyHydroGeophysX.solvers.symmetrize(A: Any) → Any[source]#

Return (A + A.T) / 2.

A covariance or resolution matrix is symmetric in exact arithmetic. Products of the form J^T C J accumulate asymmetry at the rounding level, and np.random.Generator.multivariate_normal rejects the result outright, so the halves are averaged before the matrix leaves the function that built it.