Left vector multiply

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "2.7.0"
+version = "2.8.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

README.md

@@ -8,6 +8,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.
 
+## What's new in v2.8
+*   Left multiplying by a transpose or adjoint vector (e.g., `y'*A`)
+    produces a transpose or adjoint vector output, rather than a `LinearMap`.
+
 ## What's new in v2.7
 *   Potential reduction of memory allocations in multiplication of
     `LinearCombination`s, `BlockMap`s, and real- or complex-scaled `LinearMap`s.

src/LinearMaps.jl

@@ -47,17 +47,17 @@ Base.ndims(::LinearMap) = 2
 Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
 Base.length(A::LinearMap) = size(A)[1] * size(A)[2]
 
-# check dimension consistency for y = A*x and Y = A*X
-function check_dim_mul(y::AbstractVector, A::LinearMap, x::AbstractVector)
+# check dimension consistency for multiplication C = A*B
+function check_dim_mul(C, A, B)
     # @info "checked vector dimensions" # uncomment for testing
-    m, n = size(A)
-    (m == length(y) && n == length(x)) || throw(DimensionMismatch("mul!"))
-    return nothing
-end
-function check_dim_mul(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix)
-    # @info "checked matrix dimensions" # uncomment for testing
-    m, n = size(A)
-    (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!"))
+    mA, nA = size(A) # A always has two dimensions
+    mB, nB = size(B, 1), size(B, 2)
+    if mB != nA
+        throw(DimensionMismatch("left factor has dimensions ($mA,$nA), right factor has dimensions ($mB,$nB)"))
+    end
+    if size(C, 1) != mA || size(C, 2) != nB
+        throw(DimensionMismatch("result has dimensions $(size(C)), needs ($mA,$nB)"))
+    end
     return nothing
 end
 
@@ -118,6 +118,7 @@ A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector)  = mul!(y, A,
 At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
 Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)
 
+include("left.jl") # left multiplication by a transpose or adjoint vector
 include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
 include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I
 include("transpose.jl") # transposing linear maps

src/left.jl

@@ -0,0 +1,42 @@
+# left.jl
+# "special cases" related left vector multiplication like x = y'*A
+# The subtlety is that y' is a AdjointAbsVec or a TransposeAbsVec
+# which is a subtype of AbstractMatrix but it is really "vector like"
+# so we want handle it like a vector (Issue#99)
+# So this is an exception to the left multiplication rule by a AbstractMatrix
+# that usually makes a WrappedMap.
+# The "transpose" versions may be of dubious use, but the adjoint versions
+# are useful for ensuring that (y'A)*x ≈ y'*(A*x) are both scalars.
+
+import LinearAlgebra: AdjointAbsVec, TransposeAbsVec
+
+# x = y'*A ⇐⇒ x' = (A'*y)
+Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y'))
+Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y))
+
+# mul!(x, y', A)
+Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::AdjointAbsVec, A::LinearMap)
+    @boundscheck check_dim_mul(x, y, A)
+    @inbounds mul!(adjoint(x), A', y')
+    return adjoint(x)
+end
+
+Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::AdjointAbsVec,
+        A::LinearMap, α::Number, β::Number)
+    @boundscheck check_dim_mul(x, y, A)
+    @inbounds mul!(adjoint(x), A', y', α, β)
+    return adjoint(x)
+end
+
+Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::TransposeAbsVec, A::LinearMap)
+    @boundscheck check_dim_mul(x, y, A)
+    @inbounds mul!(transpose(x), transpose(A), transpose(y))
+    return transpose(x)
+end
+
+Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::TransposeAbsVec,
+         A::LinearMap, α::Number, β::Number)
+    @boundscheck check_dim_mul(x, y, A)
+    @inbounds mul!(transpose(x), transpose(A), transpose(y), α, β)
+    return transpose(x)
+end

test/left.jl

@@ -0,0 +1,52 @@
+using Test, LinearMaps
+using LinearAlgebra: mul!
+
+# y'*L is an exception to the left multiplication rule that makes a WrappedMap
+
+
+function left_tester(L::LinearMap{T}) where {T}
+    M, N = size(L)
+    A = Matrix(L)
+
+    x = rand(T, N)
+    y = rand(T, M)
+
+    # *
+    @test y' * A ≈ y' * L
+    @test (y' * A) * x ≈ y' * (L * x) # associative
+    @test transpose(y) * A ≈ transpose(y) * L
+
+    # mul!
+    b1 = y' * A
+    b2 = similar(b1)
+    mul!(b2, y', L) # 3-arg
+    @test b1 ≈ b2
+    mul!(b2, y', L, true, false) # 5-arg
+    @test b1 ≈ b2
+
+    b1 = transpose(y) * A
+    b2 = similar(b1)
+    mul!(b2, transpose(y), L) # 3-arg
+    @test b1 ≈ b2
+    mul!(b2, transpose(y), L, true, false) # 5-arg
+    @test b1 ≈ b2
+
+    true
+end
+
+
+@testset "left mul vec" begin
+    T = ComplexF32
+    M,N = 5,5
+    L = LinearMap{T}(cumsum, reverse ∘ cumsum ∘ reverse, N)
+
+    @test left_tester(L) # FunctionMap
+    @test left_tester(L'*L) # CompositeMap
+    @test left_tester(2L) # ScaledMap
+    @test left_tester(kron(L,L')) # KroneckerMap
+    @test left_tester(2L+3L') # LinearCombination
+    @test left_tester([L L]) # BlockMap
+
+    W = LinearMap(randn(T,5,4))
+    @test left_tester(W) # WrappedMap
+end

test/runtests.jl

@@ -28,3 +28,5 @@ include("blockmap.jl")
 include("kronecker.jl")
 
 include("conversion.jl")
+
+include("left.jl")