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
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) = \{...
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
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
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, ...
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 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
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 ...
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
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
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
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
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
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
One of the standard ways to handle singular problems like Poisson with Neumann boundary conditions which lead to a sinuglar problem $Ax=b$ is to make the obtained solution orthogonal to the kernel (which is for simple domains only the constant vector) if you use something like a Krylov subspace method.
For instance, you can solve your system with conjugated ...
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
Neumann boundary conditions will show up in the right hand side of your problem, and so you can just start your iteration with any vector you want -- it doesn't have to satisfy any particular boundary condition.
Periodic boundary conditions ultimately affect the structure of the matrix. You can choose any value for $\phi(0)$ in your initial starting guess ...
3
It is a well-known theorem that conjugate gradients is guaranteed to converge (in exact arithmetic) within $r+1$ steps for an identity-plus-rank-$r$ matrix, regardless of its size. In finite precision, you would find that in the worst case, CG would nearly stagnate for the first $r+1$ iterations, then rapidly converge towards the solution for every iteration ...
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.
3
Nobody does this. It seems like everyone who has thought about the issue has come to the conclusion that the error one incurs by not computing the residual exactly is small and that computing it every few iterations would be a waste of time.
The reason, I believe, is that even for large systems with hundreds of thousands or hundreds of millions of unknowns, ...
3
It is just a projection $1 - ZZ^T$. If you are digging around in the implementation, look at KSP_PCApplyBAorAB() and related functions.
Only top voted, non community-wiki answers of a minimum length are eligible