pypcga.PCGA#

class pypcga.PCGA(*, s_init: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]], obs: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]], cov_obs: ~covmats._covariances.CovarianceMatrix, forward_model: ~typing.Callable, Q: ~covmats._covariances.CovViaEigenFactorization, drift: ~covmats._priors.DriftMatrix | None = None, prior_s_var: float | ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] | None = None, callback: ~typing.Callable | None = None, is_line_search: bool = False, is_lm: bool = False, is_direct_solve: bool = False, random_state: int | ~numpy.random._generator.Generator | ~numpy.random.mtrand.RandomState | None = Generator(PCG64) at 0x75F28C1A7060, is_objfun_exact: bool = False, max_it_lm: int = 2, alphamax_lm: float = 1000.0, lm_smin: float | None = None, lm_smax: float | None = None, max_it_ls: int = 20, max_workers_lm: int = 2, maxiter: int = 10, ftol: float = 1e-05, ftarget: float | None = None, restol: float = 0.01, logger: ~logging.Logger | None = None, is_save_jac: bool = True, eps: float = 1e-08)[source]#

Solve inverse problem with PCGA (approx to quasi-linear method)

Every values are represented as 2D np.array.

__init__(*, s_init: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]], obs: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]], cov_obs: ~covmats._covariances.CovarianceMatrix, forward_model: ~typing.Callable, Q: ~covmats._covariances.CovViaEigenFactorization, drift: ~covmats._priors.DriftMatrix | None = None, prior_s_var: float | ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] | None = None, callback: ~typing.Callable | None = None, is_line_search: bool = False, is_lm: bool = False, is_direct_solve: bool = False, random_state: int | ~numpy.random._generator.Generator | ~numpy.random.mtrand.RandomState | None = Generator(PCG64) at 0x75F28C1A7060, is_objfun_exact: bool = False, max_it_lm: int = 2, alphamax_lm: float = 1000.0, lm_smin: float | None = None, lm_smax: float | None = None, max_it_ls: int = 20, max_workers_lm: int = 2, maxiter: int = 10, ftol: float = 1e-05, ftarget: float | None = None, restol: float = 0.01, logger: ~logging.Logger | None = None, is_save_jac: bool = True, eps: float = 1e-08) → None[source]#

Initialize the instance.

Parameters:

s_init (NDArrayFloat) – 1D array of initial control parameters, i.e., initial solution for Gauss-Newton method. In theory, the choice of s_init does not affect the estimation while total number of iterations/number of forward model runs depend on s_init. Expected shape is (\(N_{s}\)) or (\(N_{s}, 1\)).
obs (numpy.ndarray, optional) – 1D array of (noisy) measurements used for inversion.
cov_obs (covmats.CovarianceMatrix) – Covariance matrix of observed data measurement errors with dimensions (\(N_{\mathrm{obs}}\), \(N_{\mathrm{obs}}\)). Also denoted \(R\) in the literature. All covariance representations provided by covmats are supported.
forward_model (Callable) – Wrapper for forward model obs = f(s). See a template python file in each example for more information. Return shape (\(N_{\mathrm{obs}}\), \(N_{e}\)).
Q (covmats.CovViaEigenFactorization) – Eigen factorization of the Covariance matrix of the inverted parameters with shape (\(N_{s}\), \(N_{s}\)).
drift (Optional[DriftMatrix], optional) – _description_, by default None
prior_s_var (Optional[Union[float, NDArrayFloat]], optional) – _description_, by default None
callback (Optional[Callable], optional) – _description_, by default None
is_line_search (bool, optional) – Whether to use line search (add ref) if the Gauss-Newton iteration fails to lower the cost function value. It comes at the cost of extra forward model runs (not parallelized). By default False.
is_lm (bool, optional) – Whether to use Levenberg Marquard inner iterations in each Gauss-Newton iterations. It consists in inflating the covariance matrix of observed data measurement errors (cov_obs), by adding weights on its diagonal. It acts as a regularization and can help to make the objective function more convex. This comes at the price of extra system inversions and max_it_lm forward model runs. The runs can be performed in parallel since it relies on forward_model. But it is the user’s responsibility to parallelize forward_model. By default False.
is_direct_solve (bool, optional) –
Whether to solve the saddle point system (Ax = b), see eq 4.53 in (add ref), with:
- the direct approach (Cholesky factorization of A), see eq (4.57) in [Collet, 2024]
- or to use the alternative iterative Krylov subspace approach as described in :cite`saibabaEfficientMethodsLargescale2012, saibabaFastAlgorithmsGeostatistical2013`.
Using direct solve is practical if the number of observations is around 100 or less. Beyond, it is advise to rely on the The default is False (using iterative Kryloc subspace by default).
random_state (Optional[ Union[int, np.random.Generator, np.random.RandomState]]) – Pseudorandom number generator state used to generate resamples. If random_state is None (or np.random), the numpy.random.RandomState singleton is used. If random_state is an int, a new RandomState instance is used, seeded with random_state. If random_state is already a Generator or RandomState instance then that instance is used.
is_objfun_exact (bool, optional) – _description_, by default False
max_it_lm (int) – Maximum number iterations when using Levenberg Marquard regularization. Only applies if is_lm is True. By default use all available CPUs.
alphamax_lm (float, optional) – Maximum weight for LM. TODO: add the formula. By default 10.0**3.
lm_smax (Optional[float], optional) – Maximum LM solution, by default None
max_it_ls (int, optional) – Maximum number of iterations when using line search, by default 20.
maxiter (int, optional) – Maximum Gauss-Newton iterations, by default 10.
ftarget (Optional[Union[float, Callable]] = None, optional) – Target objective function (stop criterion) . The iteration stops when f^{k+1} <= fmin. If None, the stop criterion is ignored. The default is None.
ftol (float, optional) – Objective function minimum change (stop criterion). The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol. Typical values for ftol on a computer with 15 digits of accuracy in double precision are as follows: ftol = 5e-3 for low accuracy; ftol = 5e-8 for moderate accuracy; ftol = 5e-14 for extremely high accuracy. If ftol = 0, the test will stop the algorithm only if the objective function remains unchanged after one iteration. The default is 1e-5.
restol (float, optional) –
Mininmum change in the update value. Below this threshold, changes are considered unisgnificant and inversion is stopped. The change is computed as the euclidean norm of the difference between the current updated values and the values at the previous iteration scaled by the euclidean norm of the values at the previous iteration:

\[\mathrm{change} = \dfrac{\left\lVert \mathbf{s}_{\ell+1} - \mathbf{s}_{\ell} \right\rVert^{2}}{\left\lVert \mathbf{s}_{\ell} \right\rVert^{2}}\]

by default 1e-2.
logger (Optional[Logger], optional) – Logger instance. If no logger is passed, there will be no output. By default None.
is_save_jac (bool, optional) – _description_, by default False
eps (float, optional) – PCGA perturbation scalar (see eq. in …), By default 1.0e-8.

Properties

`s_init`
`obs`
`forward_model`
`Q`
`callback`
`is_line_search`
`is_lm`
`is_direct_solve`
`is_objfun_exact`
`max_it_lm`
`alphamax_lm`
`lm_smin`
`lm_smax`
`max_it_ls`
`maxiter`
`ftol`
`restol`
`ftarget`
`eps`
`max_workers`
`istate`
`random_state`
`post_diagv`
`drift`
`is_save_jac`
`cov_obs_inflation_factors`
`logger`
`HX`
`HZ`
`Hs`
`invA_as_linop`
`simul_obs_init`
`cov_obs`	Get the observation errors covariance matrix.
`d_dim`	Return the number of forecast data.
`n_internal_loops`	Return the number of internal optimization loops.
`n_obs`	Return the number of observations/forecast data.
`prior_s_var`	Get the a priori variance of the control variables.
`s_dim`	Return the length of the parameters vector.

Methods

`build_cholesky`
`build_dense_A`
`build_dense_A_from_cholesky`
`display_init_parameters`
`display_objfun`
`gauss_newton`	Gauss-newton iteration
`get_A_as_linop`
`get_cov_obs_inflation_factors`	Inflation factors used in each internal loop.
`get_dense_post_cov`	Return the dense posterior covariance matrix..
`get_eigen_post_cov`	Return the posterior covariance matrix through an Eigen factorization.
`get_invA_as_linop`	Return the low rank inverse of A as a linear operator.
`get_psi`	Get the matrix HQH^{T} + lpha R.
`get_v0`
`internal_iteration_direct`
`internal_iteration_krylov_subspace`
`is_s_violate_lm_bounds`
`jac_mat`
`jac_vect`	Jacobian times Matrix (Vectors) in Parallel perturbation interval delta determined following Brown and Saad [1990]
`line_search`
`linear_iteration`	Solve the geostatistical system using a direct solver.
`loginfo`
`objective_function`	0.5(y-h(s))^TR^{-1}(y-h(s)) + 0.5*(s-Xb)^TQ^{-1}(s-Xb)
`objective_function_ls`	simul_obs with shape (n_obs, ne)
`objective_function_no_beta`	marginalized objective w.r.t.
`objective_function_no_beta_new`	marginalized objective w.r.t.
`objective_function_reg`	s_cur with shape (N_s, Ne) beta_cur with shape (N_b, Ne) 0.5(s-Xb)^TC^{-1}(s-Xb)
`rmse`	Return the root mean square error.
`run`
`solve_cholesky`