Consider the 1D Poisson equation $$\nabla^2 u = f.$$ Using finite difference method on cell corner data and a uniform grid with ghost points, I think we can write the system of equations with Neumann BCs as: $$Au = f-Au_{BC},$$ where $$ \begin{aligned} A &= \frac{1}{\Delta x^2} \left[\begin{array}{ccccccccc} -2 & 2 & & & & & \\ 1 & -2 & 1 & & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & & 1 & -2 & 1 \\ & & & & & 2 & -2 \\ \end{array} \right] \\ Ax_{BC} &= \frac{1}{\Delta x^2} \left[\begin{array}{ccccccccc} 2 \hat{n} \theta \Delta x \\ 0 \\ \\ \vdots \\ \\ 0 \\ 2 \hat{n} \theta \Delta x \\ \end{array} \right] \end{aligned} $$

And $\hat{n},\theta$ is the outward facing normal and the prescribed slope of the derivative at the boundaries.

Notably, the system is singular, which can be addressed by removing the mean of the RHS.


I know that $A$ must be symmetric, positive definite. I've done some tests without multiplying the first and last rows by 0.5 and they seem to work fine, so my question is:

Must the first and last rows of the left and right hand side be multiplied by 0.5? In other words, are there cases where it will not work without multiplication of 0.5?


I imagine that both sides of the equation would be balanced the same with and without this multiplication.

As a final note, I know that the code / algorithm may be more clear with the multiplication, since the equation explicitly satisfies symmetry, but I'm more interested in whether the multiplication is necessary or not.

Any help is greatly appreciated.

  • $\begingroup$ Related: scicomp.stackexchange.com/q/23200/713 $\endgroup$ – Kirill Oct 24 '16 at 4:02
  • 3
    $\begingroup$ Due to the minimization properties the CG algorithm is based on, it will produce approximations to the solution of $\frac12(A+A^T)x = b$ -- which is a solution of $Ax=b$ iff $A$ is symmetric, although in your case the difference is probably very small, so you don't notice. For the algorithm to converge, of course, the symmetric part $\frac12(A+A^T)$ of $A$ needs to be positive definite. $\endgroup$ – Christian Clason Oct 24 '16 at 7:22

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 is just a clever way of computing this minimizer efficiently by reusing as much information from the previous iteration as possible.)

So, for $k=n$ (at the latest), $x^n$ minimizes the function $f(x) = \frac12(x-x^*)^TA(x-x^*)$ over $\mathbb{R}^n$, which means that the gradient has to vanish. After expanding the product and using the product rule, this yields $$0 = \nabla f(x^n) = \frac12(A+A^T)x^n - \frac12(A+A^T)x^*.\tag{2}\label{cg2}$$ If $A$ is symmetric ($A^T=A$) and $Ax^*=b$, this of course implies that $Ax^n=b$; otherwise you only get a solution to $(A+A^T)x = (A+A^T)A^{-1}b$. In your case, $A$ is nearly symmetric and fairly well-conditioned, so the difference between the solutions to these two equations is probably so small that you don't see it.

Another question is whether the numerical algorithm breaks down if $A$ is not symmetric. Looking at the algorithm, this can only happen if $x^TA x =0$ for some $x\neq 0$ (which of course is impossible for positive definite matrices by definition). Now note that you can write any matrix $A$ as $$A = \frac12 (A+A^T) + \frac12(A-A^T) =: A_{\text{sym}} + A_{\text{skew}},$$ where $A_{\text{sym}}$ is the symmetric part and $A_{\text{skew}}$ the skew-symmetric part of $A$. Furthermore, $x^TA_{\text{skew}}x=0$ and hence $x^TAx = x^TA_{\text{sym}}x$ for all $x\in\mathbb{R}^n$. So the algorithm works as long as the symmetric part $A_{\text{sym}}$ of $A$ is positive definite. Conveniently, this implies that the minimization problem \eqref{cg1} is strictly convex and hence has a unique minimizer, which is the only solution to \eqref{cg2}.

TL;DR: CG works for nonsymmetric matrices as long as the symmetric part is positive definite, and will give a reasonable approximation to the solution if the skew-symmetric part is "negligible".

  • $\begingroup$ Excellent, thank you! I imagine that numerically quantifying "negligible" is a bit arbitrary, but this still gets the main point across. Thanks again! $\endgroup$ – Charles Oct 26 '16 at 18:33
  • 1
    $\begingroup$ Glad it's useful! "Negligible" was indeed quite hand-wavy, but the error can be quantified using perturbation theorems, which relate the relative error between solutions to two systems -- as measured in an appropriate matrix norm -- to the norm of the difference and the condition number of one of the matrices. $\endgroup$ – Christian Clason Oct 28 '16 at 19:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.