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
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
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 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
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 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
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 ...
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
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
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} = \...
Only top voted, non community-wiki answers of a minimum length are eligible
