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
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 ...
13
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) = \{...
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 ...
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-...
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
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 ...
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 ...
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, ...
7
You first have to make sure that your system is solvable: this happens iff the right-hand side $b$ is orthogonal to the kernel of $A$. If $A$ has a dimension-1 kernel spanned by $v$, you need to have $v^* b=0$. If that is not the case, then go back to the modelling stage and ask yourself if what you did makes sense.
Convergence of iterative methods for ...
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 ...
6
Wolfgang Bangerth's answer already says almost everything, but another subtle detail is that GMRES/MINRES minimize the norm of the residual, i.e., $\|Ax_k-b\|$, while CG minimizes the (A-)norm of the error, i.e., $\|x_*-x_k\|_A$.
In most cases what you really want to minimize is the error, and the residual serves as an imperfect proxy: recall that by the ...
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
You've started with a singular linear system of equations $Ax=b$. As a practical matter, it's unlikely that $b$ lies exactly in the range of $A$, so at best you can find a least squares solution that minimizes
$\min \| Ax - b \|_{2}$
Because the system is singular, the null space of $A$ is non-empty, and there will be an infinite number of solutions to ...
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 ...
5
The methods are unrelated. The finite element method is a way to convert a partial differential equation into a finite dimensional problem ("discretization") so it can be solved on a computer. The Conjugate Gradient method is a method to solve linear systems such as those that arise from the finite element method.
5
In some sense it doesn't matter: In finite dimensions, all norms are equivalent, so if an algorithm conveniently has the property that proving convergence in one norm is easy, then that's what people will choose. In practice, one oftentimes wants to reduce the norm of the residual by a certain factor, and from a practical perspective, a certain reduction in ...
4
I suspect there is in general not much difference between GMRES and CG for an SPD matrix.
Let's say we are solving $ Ax = b $ with $ A $ symmetric positive definite and the starting guess $ x_0 = 0 $ and generating iterates with CG and GMRES, call them $ x_k^c $ and $ x_k^g $. Both iterative methods will be building $ x_k $ from the same Krylov space $ K_k ...
4
The problem was in the underlying coupled PDE FEM model. This model includes a potential which is not fixed in any point, so any result potential = phi + C would be a solution, where C is an integration constant. I don't know why, but the direct solver (Matlab's \ using UMFPACK) finds a solution, where the iterative solver gets in massive troubles.
I first ...
4
In short, orthogonalization of the Krylov vectors occurs with respect to the operator, but not with respect to the preconditioner.
Alright, so say we want to solve $Ax=b$ with preconditioner $B$. the preconditioned-CG iteration is basically:
\begin{align*}
\hat{v}_1=\tilde{v}_1 =& Bb\\
v_1 =& \tilde{v}_1 / c_1\\
\\
\hat{v}_i =& BAv_{i-1}\\
\...
4
It looks like this paper is combining Hessian-free with truncated Newton method.
Yes, it is.
...the approach is referred to as Hessian-free method. That is because the Hessian is never computed explicitly. Instead, the product of the Hessian and a vector is obtained using finite difference approximation.
Right. This step is analogous to the "Jacobian-...
4
Let's say the condition number of the original matrix is $k_1$, and the one after refining the mesh is $k_2=4k_1$. You then want to compare the number of iterations necessary to reach a tolerance $\varepsilon$. So, assuming you are interested in a relative tolerance, in the first case, you get
$$
\left(\frac{\sqrt{k_1}-1}{\sqrt{k_1}+1}\right)^{m_1} = \...
4
While very similar, each method is slightly different and you should definitely take this into account.
The GMRES method is the simplest, it will construct an orthogonal basis for $\mathcal{K}_k(A,b)$ and select an approximate solution that minimises the 2-norm of the residual.
The MINRES method is a variation of GMRES based on the fact that for a ...
4
The nonlinear conjugate gradient method will converge for a quadratic function in $N$ steps at most. You do not have a quadratic function, and have no guarantees on convergence. It can be helpful to reset your search direction every $N$ steps or to use a better method (i.e., Newton's method). A simple function that you could use to test would be a quadratic ...
3
Probabily the CG fails to converge because your problem is ill conditioned, the condition number of your matrix is too large. For SPD matrix A you can get the condition number calculating the eigenvalues
codition_number = largest_eigenvalue/smallest_eigenvalue
Ill conditioned means that small changes on some entries of the matrix could give large changes ...
3
This tutorial talks about recalculating the residual every 50 iterations to mitigate round-off errors.
Only top voted, non community-wiki answers of a minimum length are eligible
