Skip to content

Commit

Permalink
Make LQPackedQ's postmultiplication behavior context-dependent only w…
Browse files Browse the repository at this point in the history
…here necessary.
  • Loading branch information
Sacha0 committed Sep 26, 2017
1 parent c3a3a1c commit 1b262df
Show file tree
Hide file tree
Showing 2 changed files with 61 additions and 39 deletions.
62 changes: 42 additions & 20 deletions base/linalg/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,10 +35,10 @@ lqfact(x::Number) = lqfact(fill(x,1,1))
"""
lq(A; [thin=true]) -> L, Q
Perform an LQ factorization of `A` such that `A = L*Q`. The
default is to compute a thin factorization. The LQ factorization
is the QR factorization of `A.'`. `L` is not extended with
zeros if the explicit, square form of `Q` is requested via `thin = false`.
Perform an LQ factorization of `A` such that `A = L*Q`. The default
is to compute a "thin" (or "truncated", or "reduced") factorization. The
LQ factorization is the QR factorization of `A.'`. `L` is not extended
with zeros if the explicit, square form of `Q` is requested via `thin = false`.
"""
function lq(A::Union{Number,AbstractMatrix}; thin::Bool = true)
F = lqfact(A)
Expand Down Expand Up @@ -176,17 +176,41 @@ A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasReal} =
A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasComplex} =
LAPACK.ormlq!('R', 'C', B.factors, B.τ, A)

# out-of-place right-application of LQPackedQs
# unlike their in-place equivalents, these methods: (1) check whether the applied-to
# matrix's (A's) appropriate dimension (columns for A_*, rows for Ac_*) matches the
# number of columns (nQ) of the LQPackedQ (Q), as the underlying LAPACK routine (ormlq)
# treats the implicit Q as its (nQ-by-nQ) square form; and (2) if the preceding dimensions
# do not match, these methods check whether the appropriate dimension of A instead matches
# the number of rows of the matrix of which Q is a factor (i.e. size(Q.factors, 1)),
# and if so zero-extends A as necessary for check (1) to pass (if possible).
*(A::StridedVecOrMat, Q::LQPackedQ) = _A_mul_Bq(A_mul_B!, A, Q)
A_mul_Bc(A::StridedVecOrMat, Q::LQPackedQ) = _A_mul_Bq(A_mul_Bc!, A, Q)
function _A_mul_Bq(A_mul_Bop!::FT, A::StridedVecOrMat, Q::LQPackedQ) where FT<:Function
# out-of-place right application of LQPackedQs
#
# LQPackedQ's out-of-place multiplication behavior is context dependent. specifically,
# if the inner dimension in the multiplication is the LQPackedQ's second dimension,
# the LQPackedQ behaves like its square form. if the inner dimension in the
# multiplication is the LQPackedQ's first dimension, the LQPackedQ behaves like either
# its square form or its truncated form depending on the shape of the other object
# involved in the multiplication. we treat these cases separately.
#
# (1) the inner dimension in the multiplication is the LQPackedQ's second dimension.
# in this case, the LQPackedQ behaves like its square form.
#
function A_mul_Bc(A::StridedVecOrMat, Q::LQPackedQ)
TR = promote_type(eltype(A), eltype(Q))
return A_mul_Bc!(copy_oftype(A, TR), convert(AbstractMatrix{TR}, Q))
end
function Ac_mul_Bc(A::StridedMatrix, Q::LQPackedQ)
TR = promote_type(eltype(A), eltype(Q))
C = adjoint!(similar(A, TR, reverse(size(A))), A)
return A_mul_Bc!(C, convert(AbstractMatrix{TR}, Q))
end
#
# (2) the inner dimension in the multiplication is the LQPackedQ's first dimension.
# in this case, the LQPackedQ behaves like either its square form or its
# truncated form depending on the shape of the other object in the multiplication.
#
# these methods: (1) check whether the applied-to matrix's (A's) appropriate dimension
# (columns for A_*, rows for Ac_*) matches the number of columns (nQ) of the LQPackedQ (Q),
# and if so effectively apply Q's square form to A without additional shenanigans; and
# (2) if the preceding dimensions do not match, check whether the appropriate dimension of
# A instead matches the number of rows of the matrix of which Q is a factor (i.e.
# size(Q.factors, 1)), and if so implicitly apply Q's truncated form to A by zero extending
# A as necessary for check (1) to pass (if possible) and then applying Q's square form
#
function *(A::StridedVecOrMat, Q::LQPackedQ)
TR = promote_type(eltype(A), eltype(Q))
if size(A, 2) == size(Q.factors, 2)
C = copy_oftype(A, TR)
Expand All @@ -196,11 +220,9 @@ function _A_mul_Bq(A_mul_Bop!::FT, A::StridedVecOrMat, Q::LQPackedQ) where FT<:F
else
_rightappdimmismatch("columns")
end
return A_mul_Bop!(C, convert(AbstractMatrix{TR}, Q))
return A_mul_B!(C, convert(AbstractMatrix{TR}, Q))
end
Ac_mul_B(A::StridedMatrix, Q::LQPackedQ) = _Ac_mul_Bq(A_mul_B!, A, Q)
Ac_mul_Bc(A::StridedMatrix, Q::LQPackedQ) = _Ac_mul_Bq(A_mul_Bc!, A, Q)
function _Ac_mul_Bq(A_mul_Bop!::FT, A::StridedMatrix, Q::LQPackedQ) where FT<:Function
function Ac_mul_B(A::StridedMatrix, Q::LQPackedQ)
TR = promote_type(eltype(A), eltype(Q))
if size(A, 1) == size(Q.factors, 2)
C = adjoint!(similar(A, TR, reverse(size(A))), A)
Expand All @@ -210,7 +232,7 @@ function _Ac_mul_Bq(A_mul_Bop!::FT, A::StridedMatrix, Q::LQPackedQ) where FT<:Fu
else
_rightappdimmismatch("rows")
end
return A_mul_Bop!(C, convert(AbstractMatrix{TR}, Q))
return A_mul_B!(C, convert(AbstractMatrix{TR}, Q))
end
_rightappdimmismatch(rowsorcols) =
throw(DimensionMismatch(string("the number of $(rowsorcols) of the matrix on the left ",
Expand Down
38 changes: 19 additions & 19 deletions test/linalg/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -107,10 +107,10 @@ end
@testset "correct form of Q from lq(...) (#23729)" begin
# where the original matrix (say A) is square or has more rows than columns,
# then A's factorization's triangular factor (say L) should have the same shape
# as A independent of form (thin, square), and A's factorization's orthogonal
# factor (say Q) should be a square matrix of order of A's number of columns
# independent of form (thin, square), and L and Q should have
# multiplication-compatible shapes.
# as A independent of factorization form (truncated, square), and A's factorization's
# orthogonal factor (say Q) should be a square matrix of order of A's number of
# columns independent of factorization form (truncated, square), and L and Q
# should have multiplication-compatible shapes.
m, n = 4, 2
A = randn(m, n)
for thin in (true, false)
Expand All @@ -122,18 +122,18 @@ end
# where the original matrix has strictly fewer rows than columns ...
m, n = 2, 4
A = randn(m, n)
# ... then, for a thin factorization of A, L should be a square matrix
# ... then, for a truncated factorization of A, L should be a square matrix
# of order of A's number of rows, Q should have the same shape as A,
# and L and Q should have multiplication-compatible shapes
Lthin, Qthin = lq(A, thin = true)
@test size(Lthin) == (m, m)
@test size(Qthin) == (m, n)
@test isapprox(A, Lthin * Qthin)
# ... and, for a non-thin factorization of A, L should have the same shape as A,
# Q should be a square matrix of order of A's number of columns, and L and Q
# should have multiplication-compatible shape. but instead the L returned has
# no zero-padding on the right / is L for the thin factorization, so for
# L and Q to have multiplication-compatible shapes, L must be zero-padded
# ... and, for a square/non-truncated factorization of A, L should have the
# same shape as A, Q should be a square matrix of order of A's number of columns,
# and L and Q should have multiplication-compatible shape. but instead the L returned
# has no zero-padding on the right / is L for the truncated factorization,
# so for L and Q to have multiplication-compatible shapes, L must be zero-padded
# to have the shape of A.
Lsquare, Qsquare = lq(A, thin = false)
@test size(Lsquare) == (m, m)
Expand All @@ -148,14 +148,14 @@ end
return implicitQ, explicitQ
end

m, n = 3, 3 # thin Q 3-by-3, square Q 3-by-3
m, n = 3, 3 # truncated Q 3-by-3, square Q 3-by-3
implicitQ, explicitQ = getqs(lqfact(randn(m, n)))
@test implicitQ[1, 1] == explicitQ[1, 1]
@test implicitQ[m, 1] == explicitQ[m, 1]
@test implicitQ[1, n] == explicitQ[1, n]
@test implicitQ[m, n] == explicitQ[m, n]

m, n = 3, 4 # thin Q 3-by-4, square Q 4-by-4
m, n = 3, 4 # truncated Q 3-by-4, square Q 4-by-4
implicitQ, explicitQ = getqs(lqfact(randn(m, n)))
@test implicitQ[1, 1] == explicitQ[1, 1]
@test implicitQ[m, 1] == explicitQ[m, 1]
Expand All @@ -164,7 +164,7 @@ end
@test implicitQ[m+1, 1] == explicitQ[m+1, 1]
@test implicitQ[m+1, n] == explicitQ[m+1, n]

m, n = 4, 3 # thin Q 3-by-3, square Q 3-by-3
m, n = 4, 3 # truncated Q 3-by-3, square Q 3-by-3
implicitQ, explicitQ = getqs(lqfact(randn(m, n)))
@test implicitQ[1, 1] == explicitQ[1, 1]
@test implicitQ[n, 1] == explicitQ[n, 1]
Expand Down Expand Up @@ -200,21 +200,21 @@ end
@test Ac_mul_Bc(C, implicitQ) Ac_mul_Bc(C, explicitQ)
end
# where the LQ-factored matrix has at least as many rows m as columns n,
# Q's square and thin forms have the same shape (n-by-n). hence we expect
# Q's square and truncated forms have the same shape (n-by-n). hence we expect
# _only_ *-by-n (n-by-*) C to work in A_mul_B*(C, Q) (Ac_mul_B*(C, Q)) ops.
# and hence the n-by-n C tests above suffice.
#
# where the LQ-factored matrix has more columns n than rows m,
# Q's square form is n-by-n whereas its thin form is m-by-n.
# hence we need also test *-by-m (m-by-*) C with
# A_mul_B*(C, Q) (Ac_mul_B*(C, Q)) ops, as below via m-by-m C.
# Q's square form is n-by-n whereas its truncated form is m-by-n.
# hence we need also test *-by-m C with
# A*_mul_B(C, Q) ops, as below via m-by-m C.
mA, nA = 3, 4
implicitQ, explicitQ = getqs(lqfact(randn(mA, nA)))
C = randn(mA, mA)
zeroextCright = hcat(C, zeros(eltype(C), mA))
zeroextCdown = vcat(C, zeros(eltype(C), (1, mA)))
@test *(C, implicitQ) *(zeroextCright, explicitQ)
@test A_mul_Bc(C, implicitQ) A_mul_Bc(zeroextCright, explicitQ)
@test Ac_mul_B(C, implicitQ) Ac_mul_B(zeroextCdown, explicitQ)
@test Ac_mul_Bc(C, implicitQ) Ac_mul_Bc(zeroextCdown, explicitQ)
@test_throws DimensionMismatch A_mul_Bc(C, implicitQ)
@test_throws DimensionMismatch Ac_mul_Bc(C, implicitQ)
end

0 comments on commit 1b262df

Please sign in to comment.