From ca1b5551aee65587732e2780e27f7c5d4178f006 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 22 Mar 2019 11:08:44 +0100 Subject: [PATCH 1/7] minor fixes in multiplication with Diagonals --- stdlib/LinearAlgebra/src/diagonal.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 58126a45d7d9f..7eab220f5e1e4 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -260,20 +260,20 @@ lmul!(A::Diagonal, B::Diagonal) = Diagonal(B.diag .= A.diag .* B.diag) function lmul!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix) A = adjA.parent - return lmul!(conj(A.diag), B) + return lmul!(conj(A), B) end function lmul!(transA::Transpose{<:Any,<:Diagonal}, B::AbstractMatrix) A = transA.parent - return lmul!(A.diag, B) + return lmul!(A, B) end function rmul!(A::AbstractMatrix, adjB::Adjoint{<:Any,<:Diagonal}) B = adjB.parent - return rmul!(A, conj(B.diag)) + return rmul!(A, conj(B)) end function rmul!(A::AbstractMatrix, transB::Transpose{<:Any,<:Diagonal}) B = transB.parent - return rmul!(A, B.diag) + return rmul!(A, B) end # Get ambiguous method if try to unify AbstractVector/AbstractMatrix here using AbstractVecOrMat @@ -549,13 +549,12 @@ function svd(D::Diagonal{<:Number}) end # disambiguation methods: * of Diagonal and Adj/Trans AbsVec -*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x))) +*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, parent(x), D.diag)) *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) -*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map(*, D.diag, parent(x))) +*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, parent(x), D.diag)) *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) -# TODO: these methods will yield row matrices, rather than adjoint/transpose vectors function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true) info = 0 From a7b7a6e02b6e2938051eed61d39f1d56263b660c Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 25 Mar 2019 10:23:11 +0100 Subject: [PATCH 2/7] correct rmul!(A,D), revert changes in AdjTrans(x)*D --- stdlib/LinearAlgebra/src/diagonal.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 7eab220f5e1e4..7c48d2b3c18b5 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -172,7 +172,7 @@ end function rmul!(A::AbstractMatrix, D::Diagonal) require_one_based_indexing(A) - A .= A .* transpose(D.diag) + A .= A .* permutedims(D.diag) return A end @@ -549,10 +549,10 @@ function svd(D::Diagonal{<:Number}) end # disambiguation methods: * of Diagonal and Adj/Trans AbsVec -*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, parent(x), D.diag)) +*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t*s', D.diag, parent(x))) *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) -*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, parent(x), D.diag)) +*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> t*transpose(s), D.diag, parent(x))) *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) From 5d1bcc827be6a4122e9112c03984191903e6922b Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 25 Mar 2019 10:41:52 +0100 Subject: [PATCH 3/7] [r/l]mul!: replace conj by adjoint, add transpose --- stdlib/LinearAlgebra/src/diagonal.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 7c48d2b3c18b5..5d6a98cd043f0 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -260,20 +260,20 @@ lmul!(A::Diagonal, B::Diagonal) = Diagonal(B.diag .= A.diag .* B.diag) function lmul!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix) A = adjA.parent - return lmul!(conj(A), B) + return lmul!(adjoint(A), B) end function lmul!(transA::Transpose{<:Any,<:Diagonal}, B::AbstractMatrix) A = transA.parent - return lmul!(A, B) + return lmul!(transpose(A), B) end function rmul!(A::AbstractMatrix, adjB::Adjoint{<:Any,<:Diagonal}) B = adjB.parent - return rmul!(A, conj(B)) + return rmul!(A, adjoint(B)) end function rmul!(A::AbstractMatrix, transB::Transpose{<:Any,<:Diagonal}) B = transB.parent - return rmul!(A, B) + return rmul!(A, transpose(B)) end # Get ambiguous method if try to unify AbstractVector/AbstractMatrix here using AbstractVecOrMat From fc1f32b4d7472781ba1101dc1c39e613437a6df4 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 25 Mar 2019 10:42:34 +0100 Subject: [PATCH 4/7] add tests --- stdlib/LinearAlgebra/test/diagonal.jl | 30 ++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index 0d2123db8931e..213a0f64d6a46 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -468,10 +468,20 @@ end fullBB = copyto!(Matrix{Matrix{T}}(undef, 2, 2), BB) for (transform1, transform2) in ((identity, identity), (identity, adjoint ), (adjoint, identity ), (adjoint, adjoint ), - (identity, transpose), (transpose, identity ), (transpose, transpose) ) + (identity, transpose), (transpose, identity ), (transpose, transpose), + (identity, Adjoint ), (Adjoint, identity ), (Adjoint, Adjoint ), + (identity, Transpose), (Transpose, identity ), (Transpose, Transpose)) @test *(transform1(D), transform2(B))::typeof(D) ≈ *(transform1(Matrix(D)), transform2(Matrix(B))) atol=2 * eps() @test *(transform1(DD), transform2(BB))::typeof(DD) == *(transform1(fullDD), transform2(fullBB)) end + M = randn(T, 5, 5) + MM = [randn(T, 2, 2) for _ in 1:2, _ in 1:2] + for transform in (identity, adjoint, transpose, Adjoint, Transpose) + @test lmul!(transform(D), copy(M)) == *(transform(Matrix(D)), M) + @test rmul!(copy(M), transform(D)) == *(M, transform(Matrix(D))) + @test lmul!(transform(DD), copy(MM)) == *(transform(fullDD), MM) + @test rmul!(copy(MM), transform(DD)) == *(MM, transform(fullDD)) + end end end @@ -481,10 +491,20 @@ end end @testset "Multiplication with Adjoint and Transpose vectors (#26863)" begin - x = rand(5) - D = Diagonal(rand(5)) - @test x'*D*x == (x'*D)*x == (x'*Array(D))*x - @test Transpose(x)*D*x == (Transpose(x)*D)*x == (Transpose(x)*Array(D))*x + K = 5 + x = rand(ComplexF64, K) + D = Diagonal(rand(ComplexF64, K)) + @test x'*D == x'*Array(D) == copy(x')*D == copy(x')*Array(D) + @test x'*D*x ≈ (x'*D)*x ≈ (x'*Array(D))*x + @test transpose(x)*D*x ≈ (transpose(x)*D)*x ≈ (transpose(x)*Array(D))*x + # non-commutative eltype + x = [rand(2) for _ in 1:K] + dd = [rand(2,2) for _ in 1:K] + D = Diagonal(dd) + DM = fill(zeros(2,2), K, K); DM[diagind(DM)] .= dd + @test x'*D == x'*DM == copy(x')*D == copy(x')*DM + @test x'*D*x == (x'*D)*x == (x'*DM)*x + @test transpose(x)*D*x == (transpose(x)*D)*x == (transpose(x)*DM)*x end @testset "Triangular division by Diagonal #27989" begin From 2fcfd2ea108c350cf5851a59b2e9383c7c9fa04d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 25 Mar 2019 11:06:04 +0100 Subject: [PATCH 5/7] fix typo --- stdlib/LinearAlgebra/src/diagonal.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 5d6a98cd043f0..77c9b924cae43 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -549,10 +549,10 @@ function svd(D::Diagonal{<:Number}) end # disambiguation methods: * of Diagonal and Adj/Trans AbsVec -*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t*s', D.diag, parent(x))) +*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x))) *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) -*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> t*transpose(s), D.diag, parent(x))) +*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, D.diag, parent(x))) *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) From 562412e7fb03ca9b4fd009c6fd27575703728c3f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 2 Apr 2019 11:01:17 +0200 Subject: [PATCH 6/7] relax some tests, added more tests --- stdlib/LinearAlgebra/test/diagonal.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index 213a0f64d6a46..c5cd9cb980cbf 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -493,18 +493,22 @@ end @testset "Multiplication with Adjoint and Transpose vectors (#26863)" begin K = 5 x = rand(ComplexF64, K) + xt = transpose(x) D = Diagonal(rand(ComplexF64, K)) - @test x'*D == x'*Array(D) == copy(x')*D == copy(x')*Array(D) + @test x'*D ≈ x'*Array(D) ≈ copy(x')*D ≈ copy(x')*Array(D) + @test xt*D ≈ xt*Array(D) ≈ copy(xt)*D ≈ copy(xt)*Array(D) @test x'*D*x ≈ (x'*D)*x ≈ (x'*Array(D))*x - @test transpose(x)*D*x ≈ (transpose(x)*D)*x ≈ (transpose(x)*Array(D))*x + @test xt*D*x ≈ (xt*D)*x ≈ (xt*Array(D))*x # non-commutative eltype - x = [rand(2) for _ in 1:K] - dd = [rand(2,2) for _ in 1:K] + x = [rand(ComplexF64, 2) for _ in 1:K] + xt = transpose(x) + dd = [rand(ComplexF64, 2, 2) for _ in 1:K] D = Diagonal(dd) - DM = fill(zeros(2,2), K, K); DM[diagind(DM)] .= dd - @test x'*D == x'*DM == copy(x')*D == copy(x')*DM + DM = fill(zeros(ComplexF64, 2, 2), K, K); DM[diagind(DM)] .= dd + @test x'*D ≈ x'*DM ≈ copy(x')*D ≈ copy(x')*DM + @test xt*D ≈ xt*DM ≈ copy(xt)*D ≈ copy(xt)*DM @test x'*D*x == (x'*D)*x == (x'*DM)*x - @test transpose(x)*D*x == (transpose(x)*D)*x == (transpose(x)*DM)*x + @test xt*D*x == (xt*D)*x == (xt*DM)*x end @testset "Triangular division by Diagonal #27989" begin From e215f891df04885b67e17d1ab126ddb43a04cbe9 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 3 Apr 2019 11:02:41 +0200 Subject: [PATCH 7/7] simplify tests, strict equality --- stdlib/LinearAlgebra/test/diagonal.jl | 26 +++++++++----------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index c5cd9cb980cbf..4dc697336d75a 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -491,24 +491,16 @@ end end @testset "Multiplication with Adjoint and Transpose vectors (#26863)" begin - K = 5 - x = rand(ComplexF64, K) - xt = transpose(x) - D = Diagonal(rand(ComplexF64, K)) - @test x'*D ≈ x'*Array(D) ≈ copy(x')*D ≈ copy(x')*Array(D) - @test xt*D ≈ xt*Array(D) ≈ copy(xt)*D ≈ copy(xt)*Array(D) - @test x'*D*x ≈ (x'*D)*x ≈ (x'*Array(D))*x - @test xt*D*x ≈ (xt*D)*x ≈ (xt*Array(D))*x - # non-commutative eltype - x = [rand(ComplexF64, 2) for _ in 1:K] + x = collect(1:2) xt = transpose(x) - dd = [rand(ComplexF64, 2, 2) for _ in 1:K] - D = Diagonal(dd) - DM = fill(zeros(ComplexF64, 2, 2), K, K); DM[diagind(DM)] .= dd - @test x'*D ≈ x'*DM ≈ copy(x')*D ≈ copy(x')*DM - @test xt*D ≈ xt*DM ≈ copy(xt)*D ≈ copy(xt)*DM - @test x'*D*x == (x'*D)*x == (x'*DM)*x - @test xt*D*x == (xt*D)*x == (xt*DM)*x + A = reshape([[1 2; 3 4], zeros(Int,2,2), zeros(Int, 2, 2), [5 6; 7 8]], 2, 2) + D = Diagonal(A) + @test x'*D == x'*A == copy(x')*D == copy(x')*A + @test xt*D == xt*A == copy(xt)*D == copy(xt)*A + y = [x, x] + yt = transpose(y) + @test y'*D*y == (y'*D)*y == (y'*A)*y + @test yt*D*y == (yt*D)*y == (yt*A)*y end @testset "Triangular division by Diagonal #27989" begin