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 ...
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-...
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 ...
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, \...
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
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 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
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
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
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.
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
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
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
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
It is just a projection $1 - ZZ^T$. If you are digging around in the implementation, look at KSP_PCApplyBAorAB() and related functions.
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 ...
Only top voted, non community-wiki answers of a minimum length are eligible
