Solve the linear system
A x = b
or the least-squares problem
minimize ||Ax-b||
using Minres. This is a line-by-line translation from Matlab code available at http://www.stanford.edu/group/SOL/software/minres.htm.
K = Minres(A) ; K.solve(b)
This class implements the Minres iterative solver of Paige and Saunders. Minres solves the system of linear equations Ax = b or the least squares problem min ||Ax - b||_2^2, where A is a symmetric matrix (possibly indefinite or singular) and b is a given vector.
A may be given explicitly as a matrix or be a function such that
A(x, y)
stores in y the product Ax for any given vector x. If A is an instance of some matrix class, it should have a ‘matvec’ method such that A.matvec(x, y) has the behaviour described above.
Optional keyword arguments are:
precon optional preconditioner, given as an operator (None) shift optional shift value (0.0) show display information along the iterations (True) check perform some argument checks (True) itnlim maximum number of iterations (5n) rtol relative stopping tolerance (1.0e-12)
If precon is given, it must define a positive-definite preconditioner M = C C’. The precon operator must be such that
x = precon(y)
returns the solution x to the linear system M x = y, for any given y.
If shift != 0, minres really solves (A - shift*I)x = b or the corresponding least squares problem if shift is an eigenvalue of A.
The return values (as returned by the Matlab version of Minres) are stored in the members
x, istop, itn, rnorm, Arnorm, Anorm, Acond, ynorm
of the class, after completion of solve().
Python version: Dominique Orban, Ecole Polytechnique de Montreal, 2008, translated and adapted from the Matlab version of Minres, written by
Michael Saunders, SOL, Stanford University Sou Cheng Choi, SCCM, Stanford University
See also http://www.stanford.edu/group/SOL/software/minres.html
Given x, compute the matrix-vector product y = Ax
A pure Python/numpy implementation of the Steihaug-Toint truncated preconditioned conjugate gradient algorithm as described in
T. Steihaug, The conjugate gradient method and trust regions in large scale optimization, SIAM Journal on Numerical Analysis 20 (3), pp. 626-637, 1983.
Solve the trust-region subproblem.
Given vectors s and p and a trust-region radius radius > 0, return the positive scalar sigma such that
|| s + sigma * p || = radius
in Euclidian norm. If known, supply optional argument ss whose value should be the squared Euclidian norm of s.
Solution of a trust-region subproblem using the preconditioned Lanczos method.
Create a new instance of a PyGltrContext object, representing a context to solve the quadratic problem
min <g, d> + 1/2 <d, Hd> s.t ||d|| <= radius
where either the Hessian matrix H or a means to compute matrix-vector products with H, are to be specified later.
| Parameters : |
|
||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Keywords : |
|
See the GLTR spec sheet for more information on these parameters.
Convergence of the iteration takes place as soon as
Norm(Hd + l Md + g) <= max(Norm(g) * reltol, abstol)
Solves the quadratic trust-region subproblem whose data was specified upon initialization. During the reverse communication phase, matrix vector products with the Hessian H will be computed explicitly using the matvec method of the object H. For instance, if H is an ll_mat, or csr_mat, products will be evaluated using H.matvec(x,y).
Solves the quadratic trust-region subproblem whose data was specified upon initialization. During the reverse communication phase, matrix vector products with the Hessian H will be computed implicitly using the supplied hessprod method. Given an array v, hessprod must return an array of the same size containing the result of the multiplication H*v. For instance, if the problem is from an Ampl model called nlp, the hessprod method could be
lambda v: nlp.hprod(z, v)
for some multiplier estimates z.
A general framework for implementing projected Krylov methods. Such methods are variations on all the well-known Krylov methods to solve block augmented linear systems, i.e., linear systems of the form
[ H A^T ] [ x ] = [ c ] [ A 0 ] [ y ] [ b ],
where H and A are matrices and A^T is the transpose of A. Here, H may or may not be symmetric and may or may not have stronger properties, such as positive definiteness. It may be available as an operator only or explicitly. However, all projected Krylov methods currently require that B be available explicitly.
Such matrices arise, for example, in the solution of partial differential equations (e.g., Maxwell or Navier-Stokes) by the finite-element method. For more information on general Krylov methods and projected Krylov methods, see the references below.
This module defines the ProjectedKrylov generic class. Other modules subclass ProjectedKrylov to implement specific algorithms. Currently, the following methods are implemented
Method Module Class Projected Conjugate Gradient ppcg Ppcg Projected Bi-CGSTAB pbcgstab Pbcgstab
Other projected iterative methods may be found as part of the PyKrylov package. See http://github.com/dpo/pykrylov.
| [Kel95] | C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, Philadelphia, PA, 1995. |
| [Orb08] | D. Orban, Projected Krylov Methods for Unsymmetric Augmented Systems, Cahiers du GERAD G-2008-46, GERAD, Montreal, Canada, 2008. |
Make sure constraints are consistent and residual is satisfactory
Assemble projection matrix and factorize it
- P = [ G A^T ]
- [ A 0 ],
where G is the preconditioner, or the identity matrix if no preconditioner was given.
If rhs was specified, obtain x_feasible satisfying the constraints
This is the Solve method of the abstract projectedKrylov class. The class must be specialized and this method overridden.
An implementation of the projected conjugate gradient algorithm as described in
N.I.M. Gould, M.E. Hribar and J. Nocedal, On the Solution of Equality Constrained Quadratic Programming Problems Arising in Optimization, SIAM Journal on Scientific Computing 23 (4), pp. 1376-1395, 2001.
with the addition of an optional trust-region constraint.
Bases: nlpy.krylov.projKrylov.ProjectedKrylov
Make sure constraints are consistent and residual is satisfactory
Assemble projection matrix and factorize it
- P = [ G A^T ]
- [ A 0 ],
where G is the preconditioner, or the identity matrix if no preconditioner was given.
If rhs was specified, obtain x_feasible satisfying the constraints
If fraction-to-the-boundary rule is to be enforced, compute step length to satisfy s + t*p >= btol * cur_iter.
Given vectors s and p and a trust-region radius Delta > 0, return the positive scalar sigma such that
|| s + sigma * p || = Delta
in Euclidian norm. If known, supply optional argument ss whose value should be the squared Euclidian norm of argument s.
An implementation of the projected Bi-CGSTAB algorithm as described in
D. Orban Projected Krylov Methods for Unsymmetric Augmented Systems, Cahiers du GERAD G-2008-46, GERAD, Montreal, Canada, 2008.
Bases: projKrylov.ProjectedKrylov
Make sure constraints are consistent and residual is satisfactory
Assemble projection matrix and factorize it
- P = [ G A^T ]
- [ A 0 ],
where G is the preconditioner, or the identity matrix if no preconditioner was given.
If rhs was specified, obtain x_feasible satisfying the constraints