Skip to content

Commit 38b1991

Browse files
authored
Added mul!'s, convert's, tests; expanded README (#27)
* added mul!(::AbstractMatrix,::LinearMap,::AbstractMatrix) * added 5-argument mul! * added convert functions * added tests for that * polished README
1 parent 2b20acb commit 38b1991

File tree

3 files changed

+176
-77
lines changed

3 files changed

+176
-77
lines changed

README.md

Lines changed: 50 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,22 @@
11
# LinearMaps
22

3+
[![LinearMaps](http://pkg.julialang.org/badges/LinearMaps_0.6.svg)](http://pkg.julialang.org/?pkg=LinearMaps)
34
[![LinearMaps](http://pkg.julialang.org/badges/LinearMaps_0.7.svg)](http://pkg.julialang.org/?pkg=LinearMaps)
45
[![Build Status](https://travis-ci.org/Jutho/LinearMaps.jl.svg?branch=master)](https://travis-ci.org/Jutho/LinearMaps.jl)
56
[![License](http://img.shields.io/badge/license-MIT-brightgreen.svg?style=flat)](LICENSE.md)
67
[![Coverage Status](https://coveralls.io/repos/github/Jutho/LinearMaps.jl/badge.svg?branch=master)](https://coveralls.io/github/Jutho/LinearMaps.jl?branch=master)
78
[![codecov.io](http://codecov.io/github/Jutho/LinearMaps.jl/coverage.svg?branch=master)](http://codecov.io/github/Jutho/LinearMaps.jl?branch=master)
8-
A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently.
9-
10-
## What's new.
11-
* Fully julia v0.7 compatible; dropped compatibility for previous versions of Julia.
129

13-
* `LinearMap` is now the name of the abstract type on top (used to be `AbstractLinearMap` because you could not have a function/constructor with the same name as an abstract type in older julia versions)
10+
A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently.
1411

15-
* Specifying the `eltype` of a function to be used as linear map should now use the constructor `LinearMap{T}(f, ...)`.
12+
## What's new in v2.1.0
13+
* Fully Julia v0.7 compatible; dropped compatibility for previous versions of Julia from LinearMaps.jl v2.0.0 on.
14+
* A 5-argument version for `mul!(y, A::LinearMap, x, α=1, β=0)`, which computes `y := α * A * x + β * y` and implements the usual 3-argument `mul!(y, A, x)` for the default `α` and `β`.
15+
* Synonymous `convert(Matrix, A::LinearMap)` and `convert(Array, A::LinearMap)` functions, that call the `Matrix` constructor and return the matrix representation of `A`.
16+
* Multiplication with matrices, interpreted as a block row vector of vectors:
17+
* `mul!(Y::AbstractArray, A::LinearMap, X::AbstractArray, α=1, β=0)`: applies `A` to each column of `X` and stores the result in-place in the corresponding column of `Y`;
18+
* for the out-of-place multiplication, the approach is to compute `convert(Matrix, A * X)`; this is equivalent to applying `A` to each column of `X`. In generic code which handles both `A::AbstractMatrix` and `A::LinearMap`, the additional call to `convert` is a noop when `A` is a matrix.
19+
* Full compatibility with [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl)'s `eigs` and `svds`; previously only `eigs` was working. For more, nicely collaborating packages see the [Example](#example) section.
1620

1721
## Installation
1822

@@ -24,9 +28,11 @@ Several iterative linear algebra methods such as linear solvers or eigensolvers
2428

2529
The LinearMaps package provides the following functionality:
2630

27-
1. An abstract type `LinearMap` that shares with the `AbstractMatrix` type that it responds to the functions `size`, `eltype`, `isreal`, `issymmetric`, `ishermitian` and `isposdef`, `transpose` and `adjoint` and multiplication with a vector using both `*` or the in-place version `mul!`. Linear algebra functions that uses duck-typing for its arguments can handle `LinearMap` objects similar to `AbstractMatrix` objects, provided that they can be written using the above methods. Unlike `AbstractMatrix` types, `LinearMap` objects cannot be indexed, neither using `getindex` or `setindex!`. The type `LinearMap` acts as a general purpose constructor and allows to construct linear map objects from functions, or to wrap objects of type `AbstractMatrix` or `LinearMap`. The latter functionality is useful to (re)define the properties (`isreal`, `issymmetric`, `ishermitian`, `isposdef`) of the existing matrix or linear map.
31+
1. A `LinearMap` type that shares with the `AbstractMatrix` type that it responds to the functions `size`, `eltype`, `isreal`, `issymmetric`, `ishermitian` and `isposdef`, `transpose` and `adjoint` and multiplication with a vector using both `*` or the in-place version `mul!`. Linear algebra functions that use duck-typing for their arguments can handle `LinearMap` objects similar to `AbstractMatrix` objects, provided that they can be written using the above methods. Unlike `AbstractMatrix` types, `LinearMap` objects cannot be indexed, neither using `getindex` or `setindex!`.
2832

29-
2. A framework for combining objects of type `LinearMap` and of type `AbstractMatrix` using linear combinations, transposition and composition, where the linear map resulting from these operations is never explicitly evaluated but only its matrix vector product is defined (i.e. lazy evaluation). The matrix vector product is written to minimize memory allocation by using a minimal number of temporary vectors. There is full support for the in-place version `mul!`, which should be preferred for higher efficiency in critical algorithms. In addition, it tries to recognize the properties of combinations of linear maps. In particular, compositions such as `A'*A` for arbitrary `A` or even `A'*B*C*B'*A` with arbitrary `A` and `B` and positive definite `C` are recognized as being positive definite and hermitian. In case a certain property of the resulting `LinearMap` object is not correctly inferred, the `LinearMap` method can be called to redefine the properties.
33+
2. A single method `LinearMap` function that acts as a general purpose constructor (though it is only an abstract type) and allows to construct linear map objects from functions, or to wrap objects of type `AbstractMatrix` or `LinearMap`. The latter functionality is useful to (re)define the properties (`isreal`, `issymmetric`, `ishermitian`, `isposdef`) of the existing matrix or linear map.
34+
35+
3. A framework for combining objects of type `LinearMap` and of type `AbstractMatrix` using linear combinations, transposition and composition, where the linear map resulting from these operations is never explicitly evaluated but only its matrix vector product is defined (i.e. lazy evaluation). The matrix vector product is written to minimize memory allocation by using a minimal number of temporary vectors. There is full support for the in-place version `mul!`, which should be preferred for higher efficiency in critical algorithms. In addition, it tries to recognize the properties of combinations of linear maps. In particular, compositions such as `A'*A` for arbitrary `A` or even `A'*B*C*B'*A` with arbitrary `A` and `B` and positive definite `C` are recognized as being positive definite and hermitian. In case a certain property of the resulting `LinearMap` object is not correctly inferred, the `LinearMap` method can be called to redefine the properties.
3036

3137
## Methods
3238

@@ -55,20 +61,27 @@ The LinearMaps package provides the following functionality:
5561
5662
In summary, the keyword arguments and their default values are:
5763
58-
* `ismutating`: `false` if the function `f` accepts a single vector argument corresponding to the input, and `true` if they accept two vector arguments where the first will be mutated so as to contain the result. In both cases, the resulting `A::FunctionMap` will support both the mutating as non-mutating matrix vector multiplication. Default value is guessed based on the number of arguments for the first method in the method list of `f`; it is not possible to use `f` and `fc` where only one of the two is mutating and the other is not.
59-
* `issymmetric [=false]`: whether the function represents the multiplication with a symmetric matrix. If `true`, this will automatically enable `A'*x` and `A.'*x`.
60-
* `ishermitian [=T<:Real && issymmetric]`: whether the function represents the multiplication with a hermitian matrix. If `true`, this will automatically enable `A'*x` and `A.'*x`.
64+
* `ismutating`: `false` if the function `f` accepts a single vector argument corresponding to the input, and `true` if it accepts two vector arguments where the first will be mutated so as to contain the result. In both cases, the resulting `A::FunctionMap` will support both the mutating and non-mutating matrix vector multiplication. Default value is guessed based on the number of arguments for the first method in the method list of `f`; it is not possible to use `f` and `fc` where only one of the two is mutating and the other is not.
65+
* `issymmetric [=false]`: whether the function represents the multiplication with a symmetric matrix. If `true`, this will automatically enable `A' * x` and `transpose(A) * x`.
66+
* `ishermitian [=T<:Real && issymmetric]`: whether the function represents the multiplication with a hermitian matrix. If `true`, this will automatically enable `A' * x` and `transpose(A) * x`.
6167
* `isposdef [=false]`: whether the function represents the multiplication with a positive definite matrix.
6268
63-
* `Base.Array(linearmap)`
69+
* `Base.Array(A::LinearMap)`, `Base.Matrix(A::LinearMap)`, `Base.convert(Matrix, A::LinearMap)` and `Base.convert(Array, A::LinearMap)`
70+
71+
Create a dense matrix representation of the `LinearMap` object, by multiplying it with the successive basis vectors. This is mostly for testing purposes or if you want to have the explicit matrix representation of a linear map for which you only have a function definition (e.g. to be able to use its `transpose` or `adjoint`). This way, one may conveniently make `A` act on the columns of a matrix `X`, instead of interpreting `A * X` as a composed linear map: `Matrix(A * X)`.
6472
65-
Creates a dense matrix representation of the `linearmap` object, by multiplying it with the successive basis vectors. This is mostly for testing purposes or if you want to have the explicit matrix representation of a linear map for which you only have a function definition (e.g. to be able to use its `transpose` or `adjoint`).
73+
* `SparseArrays.sparse(::LinearMap)`
6674
67-
* `SparseArrays.sparse(linearmap)`
75+
Create a sparse matrix representation of the `LinearMap` object, by multiplying it with the successive basis vectors. This is mostly for testing purposes or if you want to have the explicit sparse matrix representation of a linear map for which you only have a function definition (e.g. to be able to use its `transpose` or `adjoint`).
6876
69-
Creates a sparse matrix representation of the `linearmap` object, by multiplying it with the successive basis vectors. This is mostly for testing purposes or if you want to have the explicit sparse matrix representation of a linear map for which you only have a function definition (e.g. to be able to use its `transpose` or `adjoint`).
77+
* The following matrix multiplication methods.
7078
71-
* All matrix multiplication methods and the corresponding mutating versions.
79+
* `A * x`: applies `A` to `x` and returns the result;
80+
* `mul!(y::AbstractVector, A::LinearMap, x::AbstractVector)`: applies `A` to `x` and stores the result in `y`;
81+
* `mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix)`: applies `A` to each column of `X` and stores the results in the corresponding columns of `Y`;
82+
* `mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=1, β::Number=0)`: computes `α * A * x + β * y` and stores the result in `y`. Analogously for `X,Y::AbstractMatrix`.
83+
84+
Applying the adjoint or transpose of `A` (if defined) to `x` works exactly as in the usual matrix case: `transpose(A) * x` and `mul!(y, A', x)`, for instance.
7285
7386
## Types
7487
@@ -94,31 +107,41 @@ None of the types below need to be constructed directly; they arise from perform
94107
95108
Used to add and multiply `LinearMap` objects, don't need to be constructed explicitly.
96109
97-
## Examples
110+
## Example
98111
99-
The `LinearMap` object combines well with the iterative eigensolver `eigs` from [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl) and with iterative solvers from [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl).
112+
The `LinearMap` object combines well with methods of the following packages:
113+
* [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl): iterative eigensolver `eigs` and SVD `svds`;
114+
* [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl): iterative solvers, eigensolvers, and SVD;
115+
* [TSVD.jl](https://github.com/andreasnoack/TSVD.jl): truncated SVD `tsvd`.
100116
101117
```
102118
using LinearMaps
119+
import Arpack, IterativeSolvers, TSVD
103120

104121
function leftdiff!(y::AbstractVector, x::AbstractVector) # left difference assuming periodic boundary conditions
105-
N=length(x)
106-
length(y)==N || throw(DimensionMismatch())
107-
@inbounds for i=1:N
108-
y[i]=x[i]-x[mod1(i-1,N)]
122+
N = length(x)
123+
length(y) == N || throw(DimensionMismatch())
124+
@inbounds for i=1:N
125+
y[i] = x[i] - x[mod1(i-1, N)]
109126
end
110127
return y
111128
end
112129

113130
function mrightdiff!(y::AbstractVector, x::AbstractVector) # minus right difference
114-
N=length(x)
115-
length(y)==N || throw(DimensionMismatch())
131+
N = length(x)
132+
length(y) == N || throw(DimensionMismatch())
116133
@inbounds for i=1:N
117-
y[i]=x[i]-x[mod1(i+1,N)]
134+
y[i] = x[i] - x[mod1(i+1, N)]
118135
end
119136
return y
120137
end
121138

122-
D=LinearMap(leftdiff!, mrightdiff!, 100; ismutating=true)
123-
eigs(D'*D;nev=3,which=:SR)
139+
D = LinearMap(leftdiff!, mrightdiff!, 100; ismutating=true) # by default has eltype(D) = Float64
140+
141+
Arpack.eigs(D'D; nev=3, which=:SR) # note that D'D is recognized as symmetric => real eigenfact
142+
Arpack.svds(D; nsv=3)
143+
144+
Σ, L = IterativeSolvers.svdl(D; nsv=3)
145+
146+
TSVD.tsvd(D, 3)
124147
```

src/LinearMaps.jl

Lines changed: 34 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -21,12 +21,40 @@ Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objec
2121
Base.length(A::LinearMap) = size(A)[1] * size(A)[2]
2222

2323
Base.:(*)(A::LinearMap, x::AbstractVector) = mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
24-
function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector)
24+
function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap{T}, x::AbstractVector, α::Number=one(T), β::Number=zero(T)) where {T}
2525
length(y) == size(A, 1) || throw(DimensionMismatch("mul!"))
26-
A_mul_B!(y, A, x)
26+
if α == one(α)
27+
β == zero(β) && (A_mul_B!(y, A, x); return y)
28+
β == one(β) && (y .+= A * x; return y)
29+
# β != 0, 1
30+
rmul!(y, β)
31+
y .+= A * x
32+
return y
33+
elseif α == zero(α)
34+
β == zero(β) && (fill!(y, zero(eltype(y))); return y)
35+
β == one(β) && return y
36+
# β != 0, 1
37+
rmul!(y, β)
38+
return y
39+
else # α != 0, 1
40+
β == zero(β) && (A_mul_B!(y, A, x); rmul!(y, α); return y)
41+
β == one(β) && (y .+= rmul!(A * x, α); return y)
42+
# β != 0, 1
43+
rmul!(y, β)
44+
y .+= rmul!(A * x, α)
45+
return y
46+
end
47+
end
48+
# the following is of interest in, e.g., subspace-iteration methods
49+
function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap{T}, X::AbstractMatrix, α::Number=one(T), β::Number=zero(T)) where {T}
50+
(size(Y, 1) == size(A, 1) && size(X, 1) == size(A, 2) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!"))
51+
@inbounds @views for i = 1:size(X, 2)
52+
mul!(Y[:, i], A, X[:, i], α, β)
53+
end
54+
return Y
2755
end
2856

29-
A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
57+
A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
3058
At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
3159
Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)
3260

@@ -36,7 +64,7 @@ function Base.Matrix(A::LinearMap)
3664
T = eltype(A)
3765
mat = Matrix{T}(undef, (M, N))
3866
v = fill(zero(T), N)
39-
for i = 1:N
67+
@inbounds for i = 1:N
4068
v[i] = one(T)
4169
mul!(view(mat, :, i), A, v)
4270
v[i] = zero(T)
@@ -45,6 +73,8 @@ function Base.Matrix(A::LinearMap)
4573
end
4674

4775
Base.Array(A::LinearMap) = Base.Matrix(A)
76+
Base.convert(::Type{Matrix}, A:: LinearMap) = Matrix(A)
77+
Base.convert(::Type{Array}, A:: LinearMap) = Matrix(A)
4878

4979
# sparse: create sparse matrix representation of LinearMap
5080
function SparseArrays.sparse(A::LinearMap)

0 commit comments

Comments
 (0)