22
The associated cost of BFGS may be brought more in line with CG if you use the limited memory variants rather than the full-storage BFGS. This computes the BFGS update for the last $m$ updates efficiently by a series of rank-one updates without needing to store more than the last $m$ solutions and gradients.
In my experience, BFGS with a lot of updates ...
21
In general, all Krylov methods essentially seek a polynomial that is small when evaluated on the spectrum of the matrix. In particular, the $n$th residual of a Krylov method (with zero initial guess) can be written in the form
$$ r_n = P_n (A) b $$
where $P_n$ is some monic polynomial of degree $n$ .
If $A$ is diagonalizable, with $A=V\Lambda V^{-1}$, we ...
18
On norms
As an addendum to Reid.Atcheson's answer, I would like to clarify some issues regarding norms. At the $n^{\mathrm{th}}$ iteration, GMRES finds the polynomial $P_n$ that minimizes the $2$-norm of the residual
$$r_n = A x_n - b = \big(P_n(A) - 1 \big)b - b = P_n(A) b .$$
Suppose $A$ is SPD, so $A$ induces a norm and so does $A^{-1}$. Then
$$\begin{...
15
The matrix you're referring to is positive definite.
The eigenvalues of the matrix must be real, because symetric matrices are equal to their own conjugate transpose, and are thus Hermitian. All eigenvalues of Hermitian matrices are real.
If all entries of the matrix are positive or zero and the matrix is weakly diagonally dominant, then all eigenvalues of ...
15
In low dimensions, a well implemented BFGS method is generally both faster and more robust than CG, especially if the function is not very far from a quadratic.
Neither BFGS nor CG need any assumption about convexity; only the initial Hessian approximation (for BFGS) resp. the preconditioner (for CG) must be positive definite. But these can always be ...
15
I'll begin with a general remark: first-order information (i.e., using only gradients, which encode slope) can only give you directional information: It can tell you that the function value decreases in the search direction, but not for how long. To decide how far to go along the search direction, you need extra information (gradient descent with constant ...
14
J.M. is right about storage. BFGS requires an approximate Hessian, but you can initialize it with the identity matrix and then just calculate the rank-two updates to the approximate Hessian as you go, as long as you have gradient information available, preferably analytically rather than through finite differences. BFGS is a quasi-Newton method, and will ...
14
Gradiant descent and the conjugate gradient method are both algorithms for minimizing nonlinear functions, that is, functions like the Rosenbrock function
$
f(x_1,x_2) = (1-x_1)^2 + 100(x_2 - x_1^2)^2
$
or a multivariate quadratic function (in this case with a symmetric quadratic term)
$
f(x) = \frac{1}{2} x^T A^T A x - b^T A x.
$
Both algorithms are ...
14
First, you should forget about solving linear equations for the moment -- that's a different context from optimization, and trying to consider both on equal footing will only lead to confusion. Why some people want to motivate an optimization algorithm by mentioning a linear solver is a question only they can answer, but I would just ignore that and focus on ...
12
Iterative methods in a nutshell:
Stationary methods are in essence fixed point iterations: To solve $Ax=b$, you pick an invertible matrix $C$ and find a fixed point of
$$
x = x + Cb- CAx
$$
This converges by Banach's fixed point theorem if $\|I-CA\|<1$. The various methods then correspond to a specific choice of $C$ (e.g., for Jacobi iteration, $C=D^{-1}...
12
The conjugate gradient method is good for finding the minimum of a strictly convex functional. This is typical when you reformulate a nonlinear elliptic PDE as an optimization problem. If you want to learn about it, I recommend you read about the CG method for linear systems first, for which An Introduction to the Conjugate Gradient Method Without the ...
12
Just to clarify notation, I'll be discussing the Gauss-Newton method for the problem
$\min \phi(x)=(1/2) \| F(x) \|_{2}^{2}$
with the search direction $p$ computed as the solution to the linear system of equations
$J(x^{(k)})^{T}J(x^{(k)}) p = - J(x^{(k)})^{T} F(x^{(k)})$
where $J(x)$ is the matrix of partial derivatives of components of $F(x)$ with ...
11
The conjugate gradient method is the provably fastest iterative solver, but only for symmetric, positive-definite systems. What would be awfully convenient is if there was an iterative method with similar properties for indefinite or non-symmetric matrices.
The CG method seeks approximate solutions at each step $k$ within the Krylov subspace
$K_k(A,b) = \{...
10
Computing the determinant of a sparse matrix is typically as expensive as a direct solve, and I am skeptical that CG would be of much help in computing it. It would be possible to run CG for $n$ iterations (where $A$ is $n \times n$) in order to generate information for the entire spectrum of $A$, and to then compute the determinant as the product of the ...
9
In this context, both methods can be thought of as minimization problems of the function:
$$
\phi(\boldsymbol{x}) = \frac{1}{2}\boldsymbol{x}^T\boldsymbol{A}\boldsymbol{x} - \boldsymbol{x}^T\boldsymbol{b}
$$
When $\boldsymbol{A}$ is symmetric, then $\phi$ is minimized when $\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}$.
Gradient descent is the method ...
9
The answer is a resounding yes. The convergence rate bound of $(\sqrt{\kappa}-1) / (\sqrt{\kappa}+1)$ is sharp over the set of symmetric positive definite matrices with condition number $\kappa$. In other words, knowing nothing more about $A$ than its condition number, CG really can take $\sim\sqrt{\kappa}$ iterations to converge. Loosely speaking, the upper-...
8
One thing that CG has in its favor is that it's not minimizing the discrete $l^2$ norm for its residual polynomials (what GMRES does). It's minimizing a matrix-induced norm instead, and very often this matrix-induced norm ends up being very close to the energy norm for discretizations of physical problems, and frequently this is a much more reasonable norm ...
8
If you introduce a variable $\mathbf z=(\mathbf x^T, \mathbf y^T)^T \in \mathbb{R}^{2n}$, then you can write $f(\mathbf x, \mathbf y)=f(\mathbf z)$ and it fits the exact format you wrote in your outline of the first method and it will all work as described. You just have to match things like $\nabla_z f(\mathbf z) = (\nabla_x f(\mathbf x,\mathbf y)^T, \...
8
Eigen 3 is a nice C++ template library some of whose routines are parallelized.
c.f. Eigen documentation
The parallelization is OMP only, so if you intend to parallelise using MPI (and OMP) it is probably not suitable for your purpose.
The nice feature of Eigen is that you can swap in a high performance BLAS library (like MKL or OpenBLAS) for some routines ...
7
If your matrix is symmetric, positive definite, the CG method may converge slowly, but it converges for $n\to\infty$. The only reason it does not converge on a computer are round-off errors, in particular if the condition number of the matrix, the quotient of largest and smallest eigenvalue is large.
Experience is, that with double precision arithmetic, ...
6
If you know the null space, you can make the right hand side compatible and have the Krylov method prevent the preconditioner from causing pollution, see Why is pinning a point to remove a null space bad? for further discussion. In PETSc, this is done using the MatNullSpace object. Note that you can provide your own function to project out the null space, ...
6
The standard way is to add the constraint $u(x_0)=0$ for an arbitrarily chosen node $x_0$. This makes sure that your body can't translate or rotate and therefore takes away the zero eigenvalue. The resulting system with this constraint is positive definite even without your penalty term.
6
Others have already noted this but I think it's still worth pointing out that the determinant is not a useful quantity almost always when your matrix is large. The problem is that large matrices are most often approximations to things that are even larger dimensional (statistical samples of large populations, finite dimensional approximations to infinite ...
6
The bound is a consequence of the following three facts:
The condition number is defined as the ratio of the largest eigenvalue $\lambda_\textrm{max}$ and the smallest eigenvalue $\lambda_\textrm{min}$.
The spectral norm $\|A\|_2$ of a symmetric positive definite matrix is $\lambda_\textrm{max}$ (since, as Arnold Neumaier pointed out, it is the square root ...
6
The approximation error is $\mathcal{O}(\varepsilon)$, where $\varepsilon$ is the parameter in the finite difference calculation. It will be good to know how this approximation affects the number of CG iterations compared to having the complete Hessian? Are there theoretical/empirical understandings about this?
Yes. In short, there are problems where CG ...
6
The CG method works by producing in every iteration a vector $x^k\in\mathbb{R}^n$ that solves the minimization problem (assuming $x^0=0$ for simplicity)
$$\min_{x\in \mathcal{K}_k} \frac12(x-x^*)^TA(x-x^*),\tag{1}\label{cg1}$$
where $x^*=A^{-1}b$ and $\mathcal{K}_k$ is a Krylov space of (apart from exceptional circumstances) dimension $k$. (The CG algorithm ...
5
accumulation point = cluster point.
If a sequence oscillates between increasingly narrow neighborhoods of two different points, such as the sequence $x_n=\frac{(-1)^nn}{n+1}$, both points are accumulation points of the sequence.
In optimization, a local solver is usually reliable in practice when it is provable that each accumulation point of the sequence ...
5
The conjugate gradient method only works to solve the system
$ A x = b $
if $A$ is symmetric and positive-definite (also works for negative definite). The reason it must be symmetric is that conjugate gradient works by minimizing (or maximizing) the function
$ f(x) = \frac{1}{2} x^T A x - b^T x $
Note that the derivative is
$ f'(x) = \frac{1}{2} A^T x ...
5
1) If you're just looking to solve the PDEs without any other optimization, then my answer would be "none of them". Algorithms that discretize partial differential equations and then solve them as algebraic equations are massively parallelizable. It is possible to solve a PDE over a billion point mesh. Algorithms for nonlinear programming have made great ...
5
I highly recommend the following read:
J.R. Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain"
In short, if the matrix is non-positive definite, there is no guarantee that CG will fail. It might be able to solve it (for some RHSs and certain tolerances), but it is just not supposed to and, probably, converges slower ...
Only top voted, non community-wiki answers of a minimum length are eligible
