Reproduce 2003 spectral amge paper example

[aprokop]

In particular, we consider the Poisson equation discretized with bilinear functions
on a 32 × 32 square grid of square elements. Homogeneous Dirichlet boundary conditions
are assumed, and they are eliminated from the matrix during discretization.
Standard 2 × 2 agglomeration is forced on each level.

Unstaggered elements

metric
number of levels	2
convergence rate	0.04
grid complexity	1.47
operator complexity	1.9

[aprokop]

Steps.

1. Reproduce grid complexity Right now we get 1.23. For 2x2 element agglomerate and a single eigenvector, we should get 1 coarse dof per agglomerate, which means that the total coarsening ratio is going to be 1/4 as the agglomerates share the boundaries. So this makes sense. How did they get 1.47 in the paper, though?
2. Reproduce operator complexity Related to 1.
3. Reproduce convergence rate
   • Gauss-Seidel may be sensitive to the ordering of the degrees of freedom. We likely have a different ordering from the paper, so we may never reproduce their convergence rate exactly.

[aprokop]

Here is a relevant paragraph from the paper:

Consider Poisson’s equation, discretized using bilinear functions on rectangular
elements. The use of 2 × 2 agglomerates on the finest level necessitates $m_{\tau} = 1$ in
order to maintain $c_{op} < 1$, but this is far too restrictive for rectangular elements with large aspect 
ratios. In such cases, the assumption, implicit in the estimate, that $c_{op}$ grows at the same rate on
all coarser grids is not generally warranted. Unfortunately, more realistic growth rate estimates for 
$c_{op}$ are difficult to derive. To avoid having separate complexity measures on different levels, we 
employ a compromise in which the number of levels in the estimate is fixed. For small p, measure 
(3.10) does not restrict $m_{\tau}$ excessively; hence, we use small p in practice. In fact, empirical 
results indicate that a 3-level cost measure (p = 3) produces the most satisfactory results.
To reflect the true cost of the algorithm without resorting to an arbitrary choice of
p, we intend to refine the operator complexity estimate in future research. We also
observe that, if we really wish to restrict $c_{op} < 1$ on all levels, larger agglomerates
could likely be used, allowing more freedom in the choice of $m_{\tau}$.

In the header to the results table, they have "tests used $\mu{op}$ with $p = 3$". I think this very well means that they used $\mu{\tau} > 1$, most likely 2. This could explain the 1.47. This also means that we may need to introduce an automatic determination of the number of eigenvectors.

[aprokop]

Unfortunately, mfmg breaks right now with 2 eigenvectors. The backtrace shows that it's inside the determination of the Trilinos sparsity pattern.

[aprokop]

We reproduced the complexities with #34, but still need to fix the convergence.

[Rombur]

Smoother	Ordering	Convergence rate	Same Eigenvectors
Gauss-Seidel	default	0.105	0.095
Symmetric Gauss-Seidel	default	0.081	0.077
Gauss-Seidel	Reverse Cuthill-McKee	0.108	0.098
Symmetric Gauss-Seidel	Reverse Cuthill-McKee	0.073	0.066
Gauss-Seidel	King	0.108	0.099
Symmetric Gauss-Seidel	King	0.068	0.061
Gauss-Seidel	Reverse minimum degree	0.0995	0.092
Symmetric Gauss-Seidel	Reverse minimum degree	0.087	0.078

Using Gauss-Seidel the convergence rate is pretty stable at 0.1 Symmetric Gauss-Seidel is more dependent on the ordering but we still can't get their result :cry:

Update: Use the same eigenvectors for all the agglomerates.

[aprokop]

Sad. Btw, did you make a mistake at minimum degree Gauss-Seidel? 0.995 does not look right.

[Rombur]

Yeah a zero is missing. Any idea what we can do to fix this?

[aprokop]

Not sure what's wrong yet. Here are some hypotheses:

Check Partition of Unity Can you check (and maybe add a test) that if you sum up all the weights from all the agglomerates, you'll get a vector of all ones? I think one issue is that using global matrix diagonal may not be fully correct as it may be affected by the boundary conditions (not sure if that happens in deal.ii), while in the paper (pg.19) they divide by the sum of agglomerate contributions and not the global diagonal matrix (they are the same for all internal nodes, but I wonder if there are variations for boundary).

Check matrix entries They say:

Homogeneous Dirichlet boundary conditions are assumed, and
they are eliminated from the matrix during discretization.

So maybe when they say 32x32 system, it was originally 34x34 and with eliminated b.c. it has the same entries on the diagonal (while I think deal.ii diagonal may have different entries on the diagonal).

Eigenvectors We know that the eigenvector corresponding to the smallest eigenvalue is $e = (1, ..., 1)^T$. However, the second eigenvalue spans a subspace (as we have a standard Poisson problem) from which we pick up an eigenvector. We may choose a different eigenvector from the paper. However, I don't see how that could affect the convergence for 2-levels (I could see it affecting 3+ levels due to smoothing on the coarse level).

I would also like to know if the eigenvectors picked for different agglomerates are the same, or if the second eigenvector vary due to some random ARPACK initializations.

RHS and LHS I don't think we know what they used as the initial guess and rhs for the run. It may be that it was randomly faster than average.

Diagonal scaling of the matrix The matrix in the paper was diagonally scaled, $\hat{A} = D^{-1/2} A D^{-1/2}$. This means that the new iteration matrix for the Gauss-Seidel, $\hat{T} = (\hat{L} + \hat{D})^{-1}\hat{A} = D^{1/2}TD^{-1/2}$. It does not affect the asymptotic convergence of Gauss-Seidel, but does produce different results after finite number of steps (which is 1 for our smoother). So it may interact slightly different with the coarse grid correction. So, I'd like for us to try to diagonally scale the original matrix to see if that affects the convergence.

[Rombur]

Check Partition of Unity Can you check (and maybe add a test) that if you sum up all the weights from all the agglomerates, you'll get a vector of all ones?

It does add to one.

Check matrix entries So maybe when they say 32x32 system, it was originally 34x34 and with eliminated b.c. it has the same entries on the diagonal (while I think deal.ii diagonal may have different entries on the diagonal).

I don't think so because they say the grid is 32x32 not the system.

Eigenvectors I would also like to know if the eigenvectors picked for different agglomerates are the same, or if the second eigenvector vary due to some random ARPACK initializations.

The second vector is different for different agglomerates

RHS and LHS I don't think we know what they used as the initial guess and rhs for the run. It may be that it was randomly faster than average.

They are using a random initial guess and the rhs is zero. It is described just below equation 4.2

Diagonal scaling of the matrix So, I'd like for us to try to diagonally scale the original matrix to see if that affects the convergence.

This is not easy to test. I can easily scale the original matrix but I don't think this is sufficient. The agglomerate matrices also need to be scaled by the same scaling factor as the global one...

[aprokop]

I don't think so because they say the grid is 32x32 not the system.

But the diagonal element they have is probably still the same for all entries. Is our deal.ii matrix the same?

The second vector is different for different agglomerates

Hmm, OK. I wonder if that is done that way, is the global vector still a good approximation? If they were independent, that would have been fine. But with shared dofs, I'm not so sure. Maybe a good idea would be to try to run an experiment where we set all the the local eigenvectors corresponding to the 2nd eigenvalue to the same vector. The intuition says that nothing should be affected.

This is not easy to test. I can easily scale the original matrix but I don't think this is sufficient. The agglomerate matrices also need to be scaled by the same scaling factor as the global one...

Hmm, you may be right. Let me think about it.

[Rombur]

Using the same eigenvectors for all the agglomerates improves the convergence rate by 0.01 So the choice of the second eigenvector definitely has an effect.

[aprokop]

Hmm, you may be right. Let me think about it.

If $Ax = \lambda x$, then $D^{-1/2}Ax = \lambda D^{-1/2} x$. Let $y = D^{-1/2}x$ ($x = D^{1/2} y$), then $D^{-1/2}AD^{1/2} y = \lambda y$. Therefore, eigenvalues are going to be the same (as matrices $A$ and $D^{-1/2}AD^{1/2}$ are spectrally equivalent), and $y$ eigenvectors are going to need to be scaled by $D^{1/2}$ globally. Therefore, there is no need to scale local matrices, just need to scale the rows of restrictor after building it the usual way. I think this is correct.

[aprokop]

Per @Rombur:

I was trying to use different geometries and I found a problem. To use arpack, we need to invert the local operator. The problem is that Laplace without Dirichlet is singular (you can add an arbitrary constant to the solution). When we use umfpack right now we are lucky and because of the round-off error we are able to invert the local matrix but if I change the geometry it doesn't work anymore. This make me think that we may get the second eigenvector wrong and this could explain the convergence. I am going to see if I can use lapack instead so I don't have to invert the matrix first but we need to understand how this is supposed to work. I saw in saamge that they use arpack so they should have had the same problem.

[aprokop]

Per @Rombur:

So if I use lapack, I can use the preconditioner on distorted meshes and on a circle. However, there was a big surprise. Like we could expect lapack and arpack give the same eigenvector for the first eigenvalue but they have different eigenvectors for the second eigenvalue (there are two eigenvectors with this eigenvalues). Now if we use one or three eigenvectors, we could the same results for both arpack and lapack. However, the convergence rate with two eigenvectors is a lot worse with lapack compared to arpack. Since we don't have the full subspace with only two eigenvectors, I guess it is a matter of luck if we have the right vectors or not. When we have three eigenvectors, we have the full subspace and we are good. So with n=2, there is always some luck involved in reproducing their results. Now, even with n=3, I can't get a convergence rate of 0.04. Depending on the parameters, I get a convergence between 0.05 and 0.08.

[aprokop]

Per @Rombur:

Hyper-cube, number of eigenvectors = 1

Smoother	Ordering	Rate (ARPACK)	Rate (LAPACK)
Gauss-Seidel	default	0.267	0.267
Symmetric Gauss-Seidel	default	0.197	0.197
Gauss-Seidel	Reverse Cuthill-McKee	0.269	0.269
Symmetric Gauss-Seidel	Reverse Cuthill-McKee	0.191	0.191
Gauss-Seidel	King	0.269	0.269
Symmetric Gauss-Seidel	King	0.191	0.191
Gauss-Seidel	Reverse minimum degree	0.262	0.262
Symmetric Gauss-Seidel	Reverse minimum degree	0.218	0.218

Hyper-cube, number of eigenvectors = 2

Smoother	Ordering	Rate (APACK)	Rate (ARPACK, same)	Rate (LAPACK)
Gauss-Seidel	default	0.105	0.095	0.213
Symmetric Gauss-Seidel	default	0.081	0.077	0.147
Gauss-Seidel	Reverse Cuthill-McKee	0.108	0.098	0.216
Symmetric Gauss-Seidel	Reverse Cuthill-McKee	0.073	0.066	0.136
Gauss-Seidel	King	0.108	0.099	0.216
Symmetric Gauss-Seidel	King	0.068	0.061	0.142
Gauss-Seidel	Reverse minimum degree	0.0995	0.092	0.208
Symmetric Gauss-Seidel	Reverse minimum degree	0.087	0.078	0.165

Hyper-cube, number of eigenvectors = 3

Smoother	Ordering	Rate (ARPACK)	Rate (LAPACK)
Gauss-Seidel	default	0.075	0.075
Symmetric Gauss-Seidel	default	0.058	0.058
Gauss-Seidel	Reverse Cuthill-McKee	0.079	0.079
Symmetric Gauss-Seidel	Reverse Cuthill-McKee	0.051	0.051
Gauss-Seidel	King	0.079	0.079
Symmetric Gauss-Seidel	King	0.051	0.051
Gauss-Seidel	Reverse minimum degree	0.071	0.071
Symmetric Gauss-Seidel	Reverse minimum degree	0.057	0.057

[aprokop]

Diagonal scaling may not help with this matrix. All rows corresponding to Dirichlet dofs have 1.333 on the diagonal, and 0 offdiagonals. All other rows have 2.666 on the diagonal, and 0 or -0.333 offdiagonal. Thus, scaling would not change anything as Dirichlet rows will converge in a single iteration, and all others will be scaled exactly the same.

[Rombur]

Let's close this. I can get a convergence rate of 0.05, instead of 0.04, by taking into account the boundary conditions in the agglomerates.