It's nice to have the user's system be a matrix so "x,y,z" makes sense as u[x,y,z]
. It is really really nice. It's easy to build and plot models which are written in this form. But sadly, tools are not all designed in this manner. The main tool I am speaking of is numerical linear algebra tools like IterativeSolvers.jl or PETSc. These are designed with a linear map with vector spaces in mind.
Let me describe some of the current workflows and it will highlight the current questions and difficulties. The standard way of doing this in other languages is what I'll call "the kron
method". Open up a book or tutorial and it's probably what they do, but @haampie shows it very clearly in this Gist:
https://gist.github.com/haampie/a5edeb5df916725687b5cfd7c3b523c7
Essentially what he does is use second_order_central_diff
with a bunch of kron
commands to build an operator A
which works with vec(u)
. Then A*u
does the differentiation u_xx + u_yy
. To translate it back to the matrix, you can think of it essentially as building the A
for reshape(A*vec(u))
. Here, A
is NM x NM
, a big sparse operator so that way algorithms designed for matrix-vector multiplication "just work". This method using kron
goes to nD using similar tricks.
Now let's take a look at the other method. In the linear operator form, you can define the linear map:
This works directly on the matrix u
. What is the size(L)
? Technically this is equivalent to above, so it's actually NM x NM
. So to make this work in a linear solver, we actually need another step:
function L(u_mat)
u = reshape(u_mat,N,M)
vec(A*u + u*A)
end
The action of this linear operator on a vector of length NM is equivalent to the one built by @haampie, but does not explicitly build the operator. Using #20, it can extend easily to nD. Also, I should note that reshape
and vec
are just views on the same piece of memory (ordered the same way), so those operations are free. In the end, this is essentially a lazy matrix-free way of writing the operator that the kron
method builds.
There are advantages and disadvantages of both routes. The first one is familiar to users. Also, if you want to interface with non-Julia tools like PETSc, you need to build the full operator and send the vec(u)
. Therefore we should definitely have a way to easily build this operator. Luckily, what the kron
method needs is exactly the 1D stencil, so building that full sparse A
from PDEOperators.jl is the right thing to do.
However, the linear map is extremely useful in some cases. For example, let's take a look at the upstream operator. We haven't built it yet, but it's the one-sided (directional) derivative with a caveat: it always have to move "upstream". So A*a(t)*u
, it's the one-sided derivative in the direction of a(t)*u
, so if the sign of that flips it should flip direction. Using the kron
method, you'd have to recompute the full sparse operator each timestep to make sure it's the right direction (or lose stability). But in the lazy operator form, we can build a A*a(t) = B(t)
"lazy upstream", and do:
function L(u_mat)
u = reshape(u_mat,N,M)
vec(A*u + u*A + B(t)*u)
end
so this is strictly more powerful when it can be used. Luckily IterativeSolvers.jl seems to allow AbstractLinearMaps
(at least in what we've tried), so we can do this.
But notice what this is missing. We can't say "what is the A*vec(u) that L(u) generates"? We don't have the tools to output that. full(L)
technically does this, but in most cases that won't actually fit into memory. But LinearMaps.jl can definitely add a sparse(L)
that builds the sparse matrix for the linear map.
So in total, after going through the utility and use cases of these tools, we are missing two essential tools:
sparse(L)
in LinearMaps.jl. It can probably be implemented similarly to full(L)
? I have no idea. But this would nicely transition the dimensional model to the sparse matrix linear operator.
- Algebraic support for the
kron
method. Currently you can write sparse(A)
for A = LinearOperator ...
. However, because of directionality, you cannot do something like sparse(Axx + Ayy)
, i.e. directly build the matrix for the 2D Laplacian. It would be nice to have an algebra that supports this somehow. One way of course is to build the L
as I described and do sparse(L)
of course, but is there something more direct?