Skip to content

Commit 75e0d03

Browse files
JuthoJutho
committed
add lq and tests
1 parent 2f9a915 commit 75e0d03

7 files changed

Lines changed: 212 additions & 18 deletions

File tree

‎src/MatrixAlgebraKit.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ export left_polar!, right_polar!
3030
export left_orth, right_orth, left_null, right_null
3131
export left_orth!, right_orth!, left_null!, right_null!
3232

33+
export Native_HouseholderQR, Native_HouseholderLQ
3334
export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert,
3435
LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations,
3536
LAPACK_DivideAndConquer, LAPACK_Jacobi

‎src/common/householder.jl‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -9,21 +9,21 @@ end
99
Base.adjoint(H::Householder) = Householder(conj(H.β), H.v, H.r)
1010

1111
function householder(x::AbstractVector, r::IndexRange = axes(x, 1), k = first(r))
12-
i = findfirst(equalto(k), r)
12+
i = findfirst(==(k), r)
1313
i == nothing && error("k = $k should be in the range r = $r")
1414
β, v, ν = _householder!(x[r], i)
1515
return Householder(β, v, r), ν
1616
end
1717
# Householder reflector h that zeros the elements A[r,col] (except for A[k,col]) upon lmul!(h,A)
1818
function householder(A::AbstractMatrix, r::IndexRange, col::Int, k = first(r))
19-
i = findfirst(equalto(k), r)
19+
i = findfirst(==(k), r)
2020
i == nothing && error("k = $k should be in the range r = $r")
2121
β, v, ν = _householder!(A[r, col], i)
2222
return Householder(β, v, r), ν
2323
end
2424
# Householder reflector that zeros the elements A[row,r] (except for A[row,k]) upon rmul!(A,h')
2525
function householder(A::AbstractMatrix, row::Int, r::IndexRange, k = first(r))
26-
i = findfirst(equalto(k), r)
26+
i = findfirst(==(k), r)
2727
i == nothing && error("k = $k should be in the range r = $r")
2828
β, v, ν = _householder!(conj!(A[row, r]), i)
2929
return Householder(β, v, r), ν

‎src/implementations/lq.jl‎

Lines changed: 82 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -270,3 +270,85 @@ function _diagonal_lq!(
270270
end
271271

272272
_diagonal_lq_null!(A::AbstractMatrix, N; positive::Bool = false) = N
273+
274+
# Native logic
275+
# -------------
276+
function lq_full!(A::AbstractMatrix, LQ, alg::Native_HouseholderLQ)
277+
check_input(lq_full!, A, LQ, alg)
278+
L, Q = LQ
279+
A === Q &&
280+
throw(ArgumentError("inplace Q not supported with native LQ implementation"))
281+
_native_lq!(A, L, Q; alg.kwargs...)
282+
return L, Q
283+
end
284+
function lq_compact!(A::AbstractMatrix, LQ, alg::Native_HouseholderLQ)
285+
check_input(lq_compact!, A, LQ, alg)
286+
L, Q = LQ
287+
A === Q &&
288+
throw(ArgumentError("inplace Q not supported with native LQ implementation"))
289+
_native_lq!(A, L, Q; alg.kwargs...)
290+
return L, Q
291+
end
292+
function lq_null!(A::AbstractMatrix, N, alg::Native_HouseholderLQ)
293+
check_input(lq_null!, A, N, alg)
294+
_native_lq_null!(A, N; alg.kwargs...)
295+
return N
296+
end
297+
298+
function _native_lq!(
299+
A::AbstractMatrix, L::AbstractMatrix, Q::AbstractMatrix;
300+
positive::Bool = true # always true regardless of setting
301+
)
302+
m, n = size(A)
303+
minmn = min(m, n)
304+
@inbounds for i in 1:minmn
305+
for j in 1:(i - 1)
306+
L[i, j] = A[i, j]
307+
end
308+
β, v, L[i, i] = _householder!(conj!(view(A, i, i:n)), 1)
309+
for j in (i + 1):size(L, 2)
310+
L[i, j] = 0
311+
end
312+
H = Householder(conj(β), v, i:n)
313+
rmul!(A, H; rows = (i + 1):m)
314+
# A[i, i] == 1; store β instead
315+
A[i, i] = β
316+
end
317+
# copy remaining rows for m > n
318+
@inbounds for j in 1:size(L, 2)
319+
for i in (minmn + 1):m
320+
L[i, j] = A[i, j]
321+
end
322+
end
323+
# build Q
324+
one!(Q)
325+
@inbounds for i in minmn:-1:1
326+
β = A[i, i]
327+
A[i, i] = 1
328+
Hᴴ = Householder(β, view(A, i, i:n), i:n)
329+
rmul!(Q, Hᴴ)
330+
end
331+
return L, Q
332+
end
333+
334+
function _native_lq_null!(A::AbstractMatrix, Nᴴ::AbstractMatrix; positive::Bool = true)
335+
m, n = size(A)
336+
minmn = min(m, n)
337+
@inbounds for i in 1:minmn
338+
β, v, ν = _householder!(conj!(view(A, i, i:n)), 1)
339+
H = Householder(conj(β), v, i:n)
340+
rmul!(A, H; rows = (i + 1):m)
341+
# A[i, i] == 1; store β instead
342+
A[i, i] = β
343+
end
344+
# build Nᴴ
345+
fill!(Nᴴ, zero(eltype(Nᴴ)))
346+
one!(view(Nᴴ, 1:(n - minmn), (minmn + 1):n))
347+
@inbounds for i in minmn:-1:1
348+
β = A[i, i]
349+
A[i, i] = 1
350+
Hᴴ = Householder(β, view(A, i, i:n), i:n)
351+
rmul!(Nᴴ, Hᴴ)
352+
end
353+
return Nᴴ
354+
end

‎src/implementations/qr.jl‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -346,7 +346,8 @@ function qr_null!(A::AbstractMatrix, N, alg::Native_HouseholderQR)
346346
end
347347

348348
function _native_qr!(
349-
A::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix
349+
A::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix;
350+
positive::Bool = true # always true regardless of setting
350351
)
351352
m, n = size(A)
352353
minmn = min(m, n)
@@ -363,14 +364,15 @@ function _native_qr!(
363364
# A[j,j] == 1; store β instead
364365
A[j, j] = β
365366
end
367+
# copy remaining columns if m < n
366368
@inbounds for j in (minmn + 1):n
367369
for i in 1:size(R, 1)
368370
R[i, j] = A[i, j]
369371
end
370372
end
371373
# build Q
372374
one!(Q)
373-
for j in minmn:-1:1
375+
@inbounds for j in minmn:-1:1
374376
β = A[j, j]
375377
A[j, j] = 1
376378
Hᴴ = Householder(conj(β), view(A, j:m, j), j:m)
@@ -379,7 +381,7 @@ function _native_qr!(
379381
return Q, R
380382
end
381383

382-
function _native_qr_null!(A::AbstractMatrix, N::AbstractMatrix)
384+
function _native_qr_null!(A::AbstractMatrix, N::AbstractMatrix; positive::Bool = true)
383385
m, n = size(A)
384386
minmn = min(m, n)
385387
@inbounds for j in 1:minmn
@@ -389,10 +391,10 @@ function _native_qr_null!(A::AbstractMatrix, N::AbstractMatrix)
389391
# A[j,j] == 1; store β instead
390392
A[j, j] = β
391393
end
392-
# build Q
394+
# build N
393395
fill!(N, zero(eltype(N)))
394396
one!(view(N, (minmn + 1):m, 1:(m - minmn)))
395-
for j in minmn:-1:1
397+
@inbounds for j in minmn:-1:1
396398
β = A[j, j]
397399
A[j, j] = 1
398400
Hᴴ = Householder(conj(β), view(A, j:m, j), j:m)

‎src/interface/decompositions.jl‎

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,14 +10,23 @@
1010
# QR, LQ, QL, RQ Decomposition
1111
# ----------------------------
1212
"""
13-
Native_HouseholderQR(; blocksize, positive = false, pivoted = false)
13+
Native_HouseholderQR()
1414
1515
Algorithm type to denote a native implementation for computing the QR decomposition of
16-
a matrix using Householder reflectors, .The diagonal elements of `R` will be non-negative
16+
a matrix using Householder reflectors. The diagonal elements of `R` will be non-negative
1717
by construction.
1818
"""
1919
@algdef Native_HouseholderQR
2020

21+
"""
22+
Native_HouseholderLQ()
23+
24+
Algorithm type to denote a native implementation for computing the LQ decomposition of
25+
a matrix using Householder reflectors. The diagonal elements of `L` will be non-negative
26+
by construction.
27+
"""
28+
@algdef Native_HouseholderLQ
29+
2130
"""
2231
LAPACK_HouseholderQR(; blocksize, positive = false, pivoted = false)
2332

‎test/lq.jl‎

Lines changed: 53 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -5,9 +5,10 @@ using StableRNGs
55
using LinearAlgebra: diag, I, Diagonal
66
using MatrixAlgebraKit: LQViaTransposedQR, LAPACK_HouseholderQR
77

8-
eltypes = (Float32, Float64, ComplexF32, ComplexF64)
8+
lapack_eltypes = (Float32, Float64, ComplexF32, ComplexF64)
9+
native_eltypes = (lapack_eltypes..., BigFloat, Complex{BigFloat})
910

10-
@testset "lq_compact! for T = $T" for T in eltypes
11+
@testset "lq_compact! for T = $T" for T in lapack_eltypes
1112
rng = StableRNG(123)
1213
m = 54
1314
for n in (37, m, 63)
@@ -114,7 +115,7 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64)
114115
end
115116
end
116117

117-
@testset "lq_full! for T = $T" for T in eltypes
118+
@testset "lq_full! for T = $T" for T in lapack_eltypes
118119
rng = StableRNG(123)
119120
m = 54
120121
for n in (37, m, 63)
@@ -208,7 +209,7 @@ end
208209
end
209210
end
210211

211-
@testset "lq_compact, lq_full and lq_null for Diagonal{$T}" for T in eltypes
212+
@testset "lq_compact, lq_full and lq_null for Diagonal{$T}" for T in native_eltypes
212213
rng = StableRNG(123)
213214
atol = eps(real(T))^(3 / 4)
214215
for m in (54, 0)
@@ -252,3 +253,51 @@ end
252253
@test N isa AbstractMatrix{T} && size(N) == (0, m)
253254
end
254255
end
256+
257+
@testset "native lq_compact! for T = $T" for T in native_eltypes
258+
rng = StableRNG(123)
259+
m = 54
260+
for n in (37, m, 63)
261+
minmn = min(m, n)
262+
A = randn(rng, T, m, n)
263+
L, Q = @constinferred lq_compact(A, Native_HouseholderLQ())
264+
@test L isa Matrix{T} && size(L) == (m, minmn)
265+
@test Q isa Matrix{T} && size(Q) == (minmn, n)
266+
@test L * Q ≈ A
267+
@test all(>=(zero(real(T))), real(diag(L)))
268+
@test isisometric(Q; side = :right)
269+
Nᴴ = @constinferred lq_null(A, Native_HouseholderLQ())
270+
@test Nᴴ isa Matrix{T} && size(Nᴴ) == (n - minmn, n)
271+
@test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3)
272+
@test isisometric(Nᴴ; side = :right)
273+
274+
Ac = similar(A)
275+
L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q), Native_HouseholderLQ())
276+
@test L2 === L
277+
@test Q2 === Q
278+
Nᴴ2 = @constinferred lq_null!(copy!(Ac, A), Nᴴ, Native_HouseholderLQ())
279+
@test Nᴴ2 === Nᴴ
280+
end
281+
end
282+
283+
@testset "lq_full! for T = $T" for T in native_eltypes
284+
rng = StableRNG(123)
285+
m = 54
286+
for n in (37, m, 63)
287+
minmn = min(m, n)
288+
A = randn(rng, T, m, n)
289+
L, Q = lq_full(A, Native_HouseholderLQ())
290+
@test L isa Matrix{T} && size(L) == (m, n)
291+
@test Q isa Matrix{T} && size(Q) == (n, n)
292+
@test L * Q ≈ A
293+
@test all(>=(zero(real(T))), real(diag(L)))
294+
@test isunitary(Q)
295+
296+
Ac = similar(A)
297+
L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q), Native_HouseholderLQ())
298+
@test L2 === L
299+
@test Q2 === Q
300+
@test L * Q ≈ A
301+
@test isunitary(Q)
302+
end
303+
end

‎test/qr.jl‎

Lines changed: 55 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,10 @@ using TestExtras
44
using StableRNGs
55
using LinearAlgebra: diag, I, Diagonal
66

7-
eltypes = (Float32, Float64, ComplexF32, ComplexF64)
7+
lapack_eltypes = (Float32, Float64, ComplexF32, ComplexF64)
8+
native_eltypes = (lapack_eltypes..., BigFloat, Complex{BigFloat})
89

9-
@testset "qr_compact! and qr_null! for T = $T" for T in eltypes
10+
@testset "qr_compact! and qr_null! for T = $T" for T in lapack_eltypes
1011
rng = StableRNG(123)
1112
m = 54
1213
for n in (37, m, 63)
@@ -99,7 +100,7 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64)
99100
end
100101
end
101102

102-
@testset "qr_full! for T = $T" for T in eltypes
103+
@testset "qr_full! for T = $T" for T in lapack_eltypes
103104
rng = StableRNG(123)
104105
m = 54
105106
for n in (37, m, 63)
@@ -176,7 +177,7 @@ end
176177
end
177178
end
178179

179-
@testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in eltypes
180+
@testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in native_eltypes
180181
rng = StableRNG(123)
181182
atol = eps(real(T))^(3 / 4)
182183
for m in (54, 0)
@@ -220,3 +221,53 @@ end
220221
@test N isa AbstractMatrix{T} && size(N) == (m, 0)
221222
end
222223
end
224+
225+
226+
@testset "native qr_compact! and qr_null! for T = $T" for T in native_eltypes
227+
rng = StableRNG(123)
228+
m = 54
229+
for n in (37, m, 63)
230+
minmn = min(m, n)
231+
A = randn(rng, T, m, n)
232+
Q, R = @constinferred qr_compact(A, Native_HouseholderQR())
233+
@test Q isa Matrix{T} && size(Q) == (m, minmn)
234+
@test R isa Matrix{T} && size(R) == (minmn, n)
235+
@test Q * R ≈ A
236+
@test all(>=(zero(real(T))), real(diag(R)))
237+
N = @constinferred qr_null(A, Native_HouseholderQR())
238+
@test N isa Matrix{T} && size(N) == (m, m - minmn)
239+
@test isisometric(Q)
240+
@test maximum(abs, A' * N) < eps(real(T))^(2 / 3)
241+
@test isisometric(N)
242+
243+
Ac = similar(A)
244+
Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R), Native_HouseholderQR())
245+
@test Q2 === Q
246+
@test R2 === R
247+
N2 = @constinferred qr_null!(copy!(Ac, A), N, Native_HouseholderQR())
248+
@test N2 === N
249+
end
250+
end
251+
252+
@testset "native qr_full! for T = $T" for T in native_eltypes
253+
rng = StableRNG(123)
254+
m = 54
255+
for n in (37, m, 63)
256+
minmn = min(m, n)
257+
A = randn(rng, T, m, n)
258+
Q, R = qr_full(A, Native_HouseholderQR())
259+
@test Q isa Matrix{T} && size(Q) == (m, m)
260+
@test R isa Matrix{T} && size(R) == (m, n)
261+
@test Q * R ≈ A
262+
@test all(>=(zero(real(T))), real(diag(R)))
263+
@test isunitary(Q)
264+
265+
Ac = similar(A)
266+
Q2 = similar(Q)
267+
Q2, R2 = @constinferred qr_full!(copy!(Ac, A), (Q, R), Native_HouseholderQR())
268+
@test Q2 === Q
269+
@test R2 === R
270+
@test Q * R ≈ A
271+
@test isunitary(Q)
272+
end
273+
end

0 commit comments

Comments
 (0)