Left vector multiply (#100)

Co-authored-by: Daniel Karrasch <daniel.karrasch@posteo.de>

Author: JeffFessler

src/LinearMaps.jl

@@ -47,17 +47,15 @@ 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)
+    (mB == nA) ||
+        throw(DimensionMismatch("left factor has dimensions ($mA,$nA), right factor has dimensions ($mB,$nB)"))
+    (size(C, 1) != mA || size(C, 2) != nB) &&
+        throw(DimensionMismatch("result has dimensions $(size(C)), needs ($mA,$nB)"))
     return nothing
 end
 
@@ -118,6 +116,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,49 @@
+# 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))
+
+# multiplication with vector/matrix
+for Atype in (AdjointAbsVec, Adjoint{<:Any,<:AbstractMatrix})
+    @eval begin
+        Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap)
+            @boundscheck check_dim_mul(x, y, A)
+            @inbounds mul!(x', A', y')
+            return x
+        end
+
+        Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap,
+                α::Number, β::Number)
+            @boundscheck check_dim_mul(x, y, A)
+            @inbounds mul!(x', A', y', conj(α), conj(β))
+            return x
+        end
+    end
+end
+for Atype in (TransposeAbsVec, Transpose{<:Any,<:AbstractMatrix})
+    @eval begin
+        Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap)
+            @boundscheck check_dim_mul(x, y, A)
+            @inbounds mul!(transpose(x), transpose(A), transpose(y))
+            return x
+        end
+
+        Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap,
+                α::Number, β::Number)
+            @boundscheck check_dim_mul(x, y, A)
+            @inbounds mul!(transpose(x), transpose(A), transpose(y), α, β)
+            return x
+        end
+    end
+end

test/left.jl

@@ -0,0 +1,65 @@
+using Test, LinearMaps, LinearAlgebra
+
+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!
+    α = rand(T); β = rand(T)
+    b1 = y' * A
+    b2 = similar(b1)
+    bt = copy(b1')'
+    # bm = Matrix(bt) # TODO: this requires a generalization of the output to AbstractVecOrMat
+    @test mul!(b2, y', L) ≈ mul!(bt, y', L)# ≈ mul!(bm, y', L)# 3-arg
+    @test mul!(b2, y', L) === b2
+    @test b1 ≈ b2 ≈ bt
+    b3 = copy(b2)
+    mul!(b3, y', L, α, β)
+    @test b3 ≈ b2*β + b1*α
+
+    b1 = transpose(y) * A
+    b2 = similar(b1)
+    bt = transpose(copy(transpose(b1)))
+    @test mul!(b2, transpose(y), L) ≈ mul!(bt, transpose(y), L)
+    @test b1 ≈ b2 ≈ bt
+    b3 = copy(b2)
+    mul!(b3, transpose(y), L, α, β)
+    @test b3 ≈ b2*β + b1*α
+
+    Y = rand(T, M, 3)
+    X = similar(Y, 3, N)
+    Xt = copy(X')'
+    @test Y'L isa LinearMap
+    @test Matrix(Y'L) ≈ Y'A
+    @test mul!(X, Y', L) ≈ mul!(Xt, Y', L) ≈ Y'A
+    @test mul!(Xt, Y', L) === Xt
+    @test mul!(copy(X), Y', L, α, β) ≈ X*β + Y'A*α
+    @test mul!(X, Y', L) === X
+    @test mul!(X, Y', L, α, β) === X
+
+    true
+end
+
+@testset "left multiplication" begin
+    T = ComplexF32
+    N = 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/linearmaps.jl

@@ -47,6 +47,15 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools
         @test @inferred mul!(copy(W), M, V) ≈ AV
         @test typeof(M * V) <: LinearMap
     end
+    
+    @testset "dimension checking" begin
+        @test_throws DimensionMismatch M * similar(v, length(v) + 1)
+        @test_throws DimensionMismatch mul!(similar(w, length(w) + 1), M, v)
+        @test_throws DimensionMismatch similar(w, length(w) + 1)' * M
+        @test_throws DimensionMismatch mul!(copy(v)', similar(w, length(w) + 1)', M)
+        @test_throws DimensionMismatch mul!(similar(W, size(W).+(0,1)), M, V)
+        @test_throws DimensionMismatch mul!(copy(W), M, similar(V, size(V).+(0,1)))
+    end
 end
 
 # new type

test/runtests.jl

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