LinearMapsAA.jl
https://github.com/JeffFessler/LinearMapsAA.jl
This package is an overlay for the package
LinearMaps.jl
that allows one to represent linear operations
(like the FFT)
as a object that appears to the user like a matrix
but internally uses user-defined fast computations
for operations, especially multiplication.
With this package,
you can write and debug code
(especially for iterative algorithms)
using a small matrix A
,
and then later replace it with a LinearMapAX
object.
The extra AA
in the package name here has two meanings.
-
LinearMapAM
is a subtype ofAbstractArray{T,2}
, i.e., conforms to the requirements of anAbstractMatrix
type. -
The package was developed in Ann Arbor, Michigan :)
An AbstractArray
must support a getindex
operation.
The maintainers of the LinearMaps.jl
package
have not wished to add getindex there,
so this package adds that feature
(without committing "type piracy").
As of v0.6
,
the package produces objects of two types:
LinearMapAM
(think "Matrix") that is a subtype ofAbstractMatrix
LinearMapAO
(think "Operator") that is not a subtype ofAbstractMatrix
- The general type
LinearMapAX
is aUnion
of both. - To convert a
LinearMapAM
to aLinearMapAO
, useredim
orLinearMapAO(A)
- To convert a
LinearMapAO
to aLinearMapAM
, useundim
.
Examples
N = 6
L = LinearMap(cumsum, y -> reverse(cumsum(reverse(y))), N)
A = LinearMapAA(L) # version with no properties
A = LinearMapAA(L, (name="cumsum",))) # version with a NamedTuple of properties
Matrix(L), Matrix(A) # both the same 6 x 6 lower triangular matrix
A.name # returns "cumsum" here
Here is a more interesting example for signal processing.
using FFTW
N = 8
A = LinearMapAA(fft, y -> N*ifft(y), (N, N), (name="fft",), T=ComplexF32)
@show A[:,2]
For more details see example/fft.jl
LinearMap
objects
Features shared with Object combinations
A LinearMapAX
object supports all of the features of a LinearMap
.
In particular, if A
and B
are both LinearMapAX
objects
of appropriate sizes,
then the following each make new LinearMapAX
objects:
-
Multiplication:
A * B
-
Linear combination:
A + B
,A - B
,3A - 7B
, -
Kronecker products:
kron(A, B)
-
Concatenation:
[A B]
[A; B]
[I A I]
[A B; 2A 3I]
etc.
Caution: currently some shorthand concatenations are unsupported,
like [I I A]
, though one can accomplish that one using
lmaa_hcat(I, I, A)
Conversions
Conversion to other data types
(may require lots of memory if A
is big):
- Convert to sparse:
sparse(A)
- Convert to dense matrix:
Matrix(A)
Avoiding memory allocations
Like LinearMap
objects,
both types of LinearMapAX
objects support mul!
for storing the results in a previously allocated output array,
to avoid new memory allocations.
The basic syntax is to replace
y = A * x
with
mul!(y, A, x)
.
To make the code look more like the math,
use the InplaceOps
package:
using InplaceOps
@! y = A * x
LinearMapsAA
Features unique to A bonus feature provided by LinearMapsAA
is that a user can include a NamedTuple
of properties
with it, and then retrieve those later
using the A.key
syntax like one would do with a struct (composite type).
The nice folks over at LinearMaps.jl
helped get me started
with this feature.
Often linear operators are associated
with some properties,
e.g.,
a wavelet transform arises
from some mother wavelet,
and it can be convenient
to carry those properties with the object itself.
Currently
the properties are lost when one combines
two or more LinearMapAA
objects by adding, multiplying, concatenating, etc.
The following features are provided
by a LinearMapAX
object
due to its getindex
support:
- Columns or rows slicing:
A[:,5]
,A[end,:]
etc. return a 1D vector - Elements:
A[4,5]
(returns a scalar) - Portions:
A[4:6,5:8]
(returns a dense matrix) - Linear indexing:
A[2:9]
(returns a 1D vector) - Convert to matrix:
A[:,:]
(if memory permits) - Convert to vector:
A[:]
(if memory permits)
Operator support
A LinearMapAO
object
represents a linear mapping
from some input array size
into some output array size
specified by the idim
and odim
options.
Here is a (simplified) example for 2D MRI,
where the operator maps a 2D input array
into a 1D output vector:
using FFTW: fft, bfft
using LinearMapsAA
embed = (v, samp) -> setindex!(fill(zero(eltype(v)),size(samp)), v, samp)
N = (128,64) # image size
samp = rand(N...) .< 0.8 # random sampling pattern
K = sum(samp) # number of k-space samples
A = LinearMapAA(x -> fft(x)[samp], y -> bfft(embed(y,samp)),
(K, prod(N)) ; prop = (name="fft",), T=ComplexF32, idim=N, odim=(K,))
x = rand(N...)
z = A' * (A * x) # result is a 2D array!
For more details see example/fft.jl
The adjoint of this LinearMapAO
object
maps a 1D vector of k-space samples
into a 2D image array.
Multiplying a M × N
matrix times a N × K
matrix
can be thought of as multiplying the matrix
by each of the K
columns,
yielding a M × K
result.
Generalizing this to higher dimensional arrays,
if A::LinearMapAO
has "input dimensions" idim=(2,3)
and "output dimensions" odim=(4,5,6)
and you do A*X
where X::AbstractArray
has dimension (2,3,7,8)
,
then the output will be an Array
of dimension (4,5,6,7,8)
.
In other words, it works block-wise.
(If you really want a new LinearMapAO
, rather than an Array
,
then you must first wrap X
in a LinearMapAO
.)
This behavior deliberately departs from the non-Matrix
like behavior
in LinearMaps
where A*X
produces a new LinearMap
.
In the spirit of such generality,
this package overloads *
for LinearAlgebra.I
(and for UniformScaling
objects more generally)
such that
I * X == X
even when X
is an array of more than two dimensions.
(The original LinearAlgebra.I
can only multiply
vectors and matrices,
which suffices for matrix algebra,
but not for general linear algebra.)
Caution:
The LinearMapAM
type should be quite stable now,
whereas LinearMapAO
is new in v0.6
.
The conversions redim
and undim
are probably not thoroughly tested.
The safe bet is to use all
LinearMapAM
objects
or all
LinearMapAO
objects
rather than trying to mix and match.
Caution
An AbstractArray
also must support a setindex!
operation
and this package provides that capability,
mainly for completeness
and as a proof of principle,
solely for the LinearMapAM
type.
Examples:
A[2,3] = 7
A[:,4] = ones(size(A,1))
A[end] = 0
A single setindex!
call is reasonably fast,
but multiple calls add layers of complexity
that are likely to slow things down.
In particular, trying to do something like the Gram-Schmidt procedure
"in place" with an AbstractArray
would be insane.
In fact, LinearAlgebra.qr!
works only with a StridedMatrix
not a general AbstractMatrix
.
Related packages
LinearOperators.jl
also provides getindex
-like features,
but slicing there always returns another operator,
unlike with a matrix.
In contrast,
a LinearMapAM
object is designed to behave
akin to a matrix,
except for operations like svd
and pinv
that are unsuitable for large-scale problems.
However, one can try
Arpack.svds(A)
to compute a few SVD components.
LazyArrays.jl
and
BlockArrays.jl
also have some related features,
but only for arrays,
not linear operators defined by functions,
so LinearMaps
is more comprehensive.
This package provides similar functionality
as the Fatrix
/ fatrix
object in the
Matlab version of MIRT.
In particular,
the odim
and idim
features of those objects
are similar to those here.
FunctionOperators.jl
supports inDims
and outDims
features.
Being a sub-type of AbstractArray
can be useful
for other purposes,
such as using the nice
Kronecker.jl
package.
Multiplication properties
It can help developers and users to know how multiplication operations should behave.
Type | Shorthand |
---|---|
LinearMapAO |
O |
LinearMapAM |
M |
LinearMap |
L |
AbstractVector |
v |
AbstractMatrix |
X |
AbstractArray |
A |
LinearAlgebra.I |
I |
For left * right
multiplication the results are as follows.
Left | Right | Result |
---|---|---|
M |
v |
v |
v' |
M |
v' |
M |
X |
X |
X |
M |
X |
M |
M |
M |
M |
L |
M |
L |
M |
M |
O |
A |
A |
A |
O |
A |
O |
O |
O |
I |
A |
A |
The following subset of the above operations also work
for the in-place version mul!(result, left, right)
:
Left | Right | Result |
---|---|---|
M |
v |
v |
v' |
M |
v' |
M |
X |
X |
X |
M |
X |
O |
A |
A |
A |
O |
A |
Credits
This software was developed at the
University of Michigan
by
Jeff Fessler
and his
group,
with substantial inspiration drawn
from the LinearMaps
package.
Compatability
- Version 0.2.0 tested with Julia 1.1 and 1.2
- Version 0.3.0 requires Julia 1.3
- Version 0.6.0 assumes Julia 1.4
- master assumes Julia 1.5
Getting started
This package is registered in the
General
registry,
so you can install it at the REPL with ] add LinearMapAA
.
Here are detailed installation instructions
This package is included in the
Michigan Image Reconstruction Toolbox
MIRT.jl
and is exported there
so that MIRT users can use it
without "separate" installation.