Solving a linear system with matrix arguments

Question

We're all familiar with the many computational methods to solve the standard linear system

$$ Ax=b. $$However, I'm curious if there are any "standard" computational methods for solving a more general (finite dimensional) linear system of the form

$$ LA=B, $$ where, say, $A$ is an $m_1\times n_1$ matrix, $B$ is an $m_2\times n_2$ matrix, and $L$ is a linear operator taking $m_1\times n_1$ matrices to $m_2\times n_2$ matrices, which does not involve vectorizing the matrices, i.e. converting everything to the standard $Ax=b$ form.

The reason I ask is I need to solve the following equation for $u$:

$$ (R^*R+\lambda I)u=f $$ where $R$ is the 2d Radon transform, $R^*$ its adjoint, and both $u$ and $f$ are 2d arrays (images). It is possible to vectorize this equation, but it's a pain, especially if we go to 3D.

More generally, what about $nD$ arrays? For instance, solving $LA=B$ where $A$ and $B$ are 3D arrays (I'll need to do this with Radon transform at some point as well).

Thanks ahead, and feel free to ship me off to another StackExchange if you feel the need.

Answer 1

Yes, you've got it right, and it will indeed work fine when you upgrade to 3-D. The easiest part, really, is the inner product---just do a standard dot product on the equivalent, unrolled $\mathbb{R}^n$ vectors. Since you'll likely have the data stored contiguously anyway, you can do this in-place. This even works with complex vector spaces---just treat the complex values as pairs of real values. That's because for CG you need the real inner product $\langle y, x\rangle\triangleq \mathop{\textrm{Re}}(y^Hx)$.

One thing you have to be careful of when you're implementing CG (or similar iterative approaches) with general linear operators is to implement the adjoint of your linear operator properly. That is, people often get $y=F(x)$ right, but make a mistake implementing $z=F^*(y)$.

I recommend implementing a simple test that takes advantage of the following identity: for any conforming $x$ and $y$, $$\langle y, F(x) \rangle = \langle F^*(y), x \rangle.$$ So what you do is generate random values of $x$ and $y$, run them through your forward and adjoint operations, respectively, and compute the two inner products above. Make sure they match to within reasonable precision, and repeat a few times.

EDIT: what do you do if your linear operator is supposed to be symmetric? Well, you do need to verify that symmetry, too. So use the same test, just noting that $F=F^*$---apply the same operation to $x$ and $y$. Of course, the OP has both an asymmetric operator and a symmetric one to deal with...

Answer 2

As it turns out, because my system is symmetric and positive definite (since my linear operator is written as $R^*R+\lambda I$), conjugate gradient can be adapted to solve this type of equation iteratively. The only modification comes when computing inner products - i.e. a typical inner product computation in CG looks like $r_k^Tr_k$ or $p_k^TAp_k$. In the modified version, we use the Frobenius inner product, which can be computed by summing the entries of the Hadamard (pointwise) product. I.e.

$$ \langle A,B\rangle=\sum_{i,j}A_{ij}B_{ij} $$

I suspect this will go through just fine when I upgrade to 3D arrays, though I have yet to see the Frobenius inner product defined on 3D arrays (I'll work under the assumption that I can again just sum the pointwise product).

I'd still be interested in more general methods if anyone knows of any!