# SPD matrices with right hand sides

I'm looking for sparse SPD matrices with right hand side? There is UF collection of sparse matrices, however, I'm not sure how do I search of the matrices of these kind efficiently (I'm doing a naive search which hasn't given me any results so far and it takes arbitrary long for some of the matrices).

• There's a search page on the website for the collection -- is this what you're using? It would help if you were more specific. – Daniel Shapero May 22 '13 at 22:30
• You mean you are looking for linear systems $Ax=b$ where $A$ is SPD and $b$ is given? What's your goal why you need such a $b$? – Wolfgang Bangerth May 22 '13 at 22:47
• If you're just interested in testing a solver, you can randomly generate an $x$ vector, multiply $Ax$ to get $b$, and then solve for $x$ using that $b$. – Brian Borchers May 23 '13 at 0:02
• @BrianBorchers Yes, but it's better to select $x$ that is not accidentally orthogonal to the smallest eigenvector of $A$. If you use a Gassian distribution with mean zero for a Laplacian, for example, the convergence with many methods will be anomalously fast (as you refine the grid) because the solution $x$ becomes more exactly orthogonal to the smallest eigenvector as the random process homogenizes. – Jed Brown May 23 '13 at 3:12
• @BrianBorchers 1. Your first comment suggested randomly generating $x$, then computing $b \gets A x$, then solving. This is convenient because you can evaluate an error, but is orthogonal to my point: 2. If your distribution is orthogonal to that lowest eigenvector, every random vector ($x$ or $b$) that you draw will have this property of converging artificially fast. People often don't remember to add correlation (shifts, rotations, etc) to the vectors they draw. This mistake gets published from time to time, so it's worth increasing awareness. – Jed Brown May 23 '13 at 14:43

We acknowledge that most matrix collections don't include physically reasonable right hand sides. If we want to test solver performance, we have to generate problems to solve --- either the right-hand sides $b$ or solutions $x$ from which we compute $b \gets A x$ and then solve $A y = b$ to some precision. The latter is convenient because it lets us check accuracy. Solvers like unpreconditioned GMRES are optimal (over a subspace) in the 2-norm of the residual, thus the $A^* A$-norm of the error. This is a significantly weaker norm than the 2-norm of the error, for example, thus it really is useful to compare solvers in the error norm.

### Dangers of choosing a random vector

A common method of generating $x$ or $b$ is to draw from an independent random distribution. This is flawed because these random vectors will be nearly orthogonal to the eigenvectors associated with the smallest eigenvalues -- exactly the modes that slow solver convergence. For example, if you use a Gaussian distribution with mean zero for a Neumann Laplacian (has a constant null space), then the average over any significant portion of the domain will be nearly zero, meaning that no long-range interactions take place. In this case, a coarse grid becomes unnecessary since you only have to move information far enough for the random process to homogenize.

If your distribution is orthogonal to that lowest eigenvector, every random vector ($x$ or $b$ doesn't matter when comparing to an eigenvector $z$, though if $z^T x$ is close to zero, $z^T b = z^T A x$ will be even smaller) that you draw will have this property of being nearly orthogonal to the "problematic modes", thus converging artificially fast. An independent random distribution will be orthogonal to most or all these "problematic modes" (e.g., translations and rotations for elasticity).

This mistake gets published from time to time, so it's worth increasing awareness.

### A Fix

If you don't know anything about the problem, it's inconvenient to tune the random distribution so that its expected value has nonzero component in the direction of all the low-energy eigenvectors---you would have to solve an expensive eigenproblem to find those "bad" vectors. Instead, I would recommend the following ad-hoc procedure:

1. Recursively bisect the graph.
2. At each level of the bisection tree, choose a random value from a random distribution with nonzero mean.
3. Define the vector $x$ to be the sum of the random value assigned to all the parts that contain it (i.e., the sum of the path from the root of the tree).

This is much like random sampling in a wavelet basis and provides correlation at all scales. It doesn't require any special knowledge about the SPD matrix other than its sparsity pattern.

### Extensions and alternatives

• Michael Grant asks "What about dense systems?" My reply: For an unstructured dense matrix, I would threshold since that should be some measure of locality. For a structured dense matrix (e.g., an $H$-matrix), you already have a hierarchy. The details of the bisection/partitioning isn't that important. All it needs to do is provide correlation at different scales so that the distribution isn't orthogonal to the problematic eigenvectors.