Skip to content

Latest commit

 

History

History
457 lines (375 loc) · 14.6 KB

README.md

File metadata and controls

457 lines (375 loc) · 14.6 KB
Raw

LinearMapsAA.jl

https://github.com/JeffFessler/LinearMapsAA.jl

action status build status pkgeval status codecov.io license docs stable docs dev

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.

As of v0.6, the package produces objects of two types:

  • LinearMapAM (think "Matrix") that is a subtype of AbstractMatrix.
  • LinearMapAO (think "Operator") that is not a subtype of AbstractMatrix.
  • The general type LinearMapAX is a Union of both.
  • To convert a LinearMapAM to a LinearMapAO, use redim or LinearMapAO(A)
  • To convert a LinearMapAO to a LinearMapAM, use undim.

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 LinearMapsAA
using FFTW: fft, bfft
N = 8
A = LinearMapAA(fft, bfft, (N, N), (name="fft",), T=ComplexF32)
@show A[:,2]

For more details see the examples in the documentation.

Features shared with LinearMap objects

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 # shorthand for mul!(y, A, x)

Features unique to LinearMapsAA

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!
typeof(A) # LinearMapAO{ComplexF32, 1, 2}

For more details see FFT example in the documentation.

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.

Here is a corresponding example (not useful; just for illustration).

using LinearMapsAA
idim = (2,3)
odim = (4,5,6)
forward = x -> repeat(reshape(x, (idim[1],1,idim[2])) ; outer=(2,5,2))
A = LinearMapAA(forward,
    (prod(odim), prod(idim)) ; prop = (name="test",), idim, odim)
x = rand(idim..., 7, 8)
y = A * x

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.

Historical note about getindex

An AbstractArray must support a getindex operation. The maintainers of the LinearMaps.jl package originally did not wish to add getindex there, so this package added that feature (while avoiding "type piracy"). Eventually, partial getindex support, specifically slicing, was added in v3.7 there. As of v0.11, this package uses that getindex implementation and also supports only slicing. This is a breaking change that could be easily reverted, so please submit an issue if you have a use case for more general use of getindex.

Historical note about setindex!

The Julia manual section on the AbstractArray interface implies that an AbstractArray should support a setindex! operation. Versions of this package prior to v0.8.0 provided that capability, mainly for completeness and as a proof of principle, solely for the LinearMapAM type. However, the reality is that many sub-types of AbstractArray in the Julia ecosystem, such as LinearAlgebra.Diagonal, understandably do not support setindex!, and it there seems to be no use for it here either. Supporting setindex! seems impossible with a concrete type for a function map, so it is no longer supported. The key code is relegated to the archive directory.

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.

LazyAlgebra.jl also has many related features, and supports nonlinear mappings.

SciML/SciMLOperators.jl seems to be designed for "matrix-free" operators that are functions of some possibly changing parameters.

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.

JOLI.jl

Being a sub-type of AbstractArray can be useful for other purposes, such as using the nice Kronecker.jl package.

Will this list keep growing, or will the community settle on some common AbstractLinearMap base?

Inter-operability

To promote inter-operability of different linear mapping packages, this package provides methods for wrapping other operator types into LinearMapAX types. The syntax is simply LinearMapAA(L; kwargs...) where L can be any of the following types currently:

Submit an issue or make a PR if there are other operator types that one would like to have supported. To minimize package dependencies, the wrapping code for a LinearOperator uses package extensions.

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

There is one more special multiplication property. If O is a LinearMapAO and xv is Vector of AbstractArrays, then O * xv is equivalent to [O * x for x in xv]. This is useful, for example, in dynamic imaging where one might store a video sequence as a vector of 2D images, rather than as a 3D array. There is also a corresponding 5-argument mul!. Each array in the Vector xv must be compatible with multiplication on the left by O.

Credits

This software was developed at the University of Michigan by Jeff Fessler and his group, with substantial inspiration drawn from the LinearMaps package.

Compatibility

  • 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
  • Version 0.7.0 assumes Julia 1.6
  • Version 0.11.0 assumes Julia 1.9

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.