Relax singularity check for tolerance

Author: Pawel Latawiec

src/lal.jl

@@ -106,7 +106,6 @@ mutable struct LookAheadLanczosDecomp{OpT, OptT, VecT, MatT, ElT, ElRT}
 
     # Estimate of norm(A), see [^Freund1993]
     nA::ElRT
-    nA_recompute::ElRT
 
     # Logs and options
     log::LookAheadLanczosDecompLog
@@ -227,7 +226,7 @@ function LookAheadLanczosDecomp(
         D, E, Flastcol, Flastrow, F̃lastcol, G, H,
         U, L,
         n, k, l, kstar, lstar, mk, nl,
-        false, false, nA, nA,
+        false, false, nA,
         ld_log,
         LookAheadLanczosDecompOptions(
             max_iter,
@@ -311,7 +310,7 @@ function _update_PQ_sequence!(ld)
     # Alg. 5.2.3
     _update_Flastrow!(ld)
     # Alg. 5.2.4
-    ld.innerp = inner_ok && _is_singular(ld.E)
+    ld.innerp = inner_ok && !isempty(ld.E) && _is_singular(last(blocks(ld.E)))
     # Alg. 5.2.5
     _update_U!(ld, ld.innerp)
     # Alg. 5.2.6
@@ -374,7 +373,7 @@ function _update_VW_sequence!(ld)
     # Alg. 5.2.16
     _update_Flastcol!(ld)
     # Alg. 5.2.17
-    ld.innerv = inner_ok && _is_singular(ld.D[ld.nl[ld.l]:end, ld.nl[ld.l]:end])
+    ld.innerv = inner_ok && _is_singular(last(blocks(ld.D)))
     # Alg. 5.2.18
     _update_L!(ld, ld.innerv)
     # Alg. 5.2.19
@@ -436,7 +435,6 @@ function _update_D!(ld)
     # Eq. 3.15, (D Γ)ᵀ = (D Γ)
     # D[n, n] = wtv
 
-    # TODO: closed block
     if isone(ld.n) || _VW_block_size(ld) == 1
         _start_new_block!(ld.D, ld.wtv)
     else
@@ -464,7 +462,6 @@ function _update_Flastrow!(ld)
     # Note: D is at D_n, L is at L_{n-1}
     # Eq. 5.2 (w/ indices advanced): 
     # F_{n} = D_{n}L[1:n, 1:n] + l[n+1, n]D_{n}[1:n, n+1][0 ... 0 1]
-    # TODO: block
     if !isone(ld.n) # We only need to do this if we are constructing a block
         ld.Flastrow = reshape(ld.D[end:end, :] * ld.L, :)
         ld.F̃lastcol = ld.Flastrow .* ld.γ[1:end-1] ./ ld.γ[end]
@@ -589,7 +586,9 @@ function _append_PQ!(ld)
 end
 
 # 5.2.4, 5.2.17
-_is_singular(A) = !isempty(A) && minimum(svdvals(A)) < eps(real(eltype(A)))
+# Freund, R. W., Gutknecht, M. H., & Nachtigal, N. M. (1993). An Implementation of the Look-Ahead Lanczos Algorithm for Non-Hermitian Matrices Part I. SIAM Journal on Scientific Computing, 14(1), 137–158. https://doi.org/10.1137/0914009
+# suggests eps()^(1/3)
+_is_singular(A) = !isempty(A) && minimum(svdvals(A)) < eps(real(eltype(A)))^(1/3)
 
 function _update_E!(ld)
     # E = QtAP

src/limited_memory_matrices.jl

@@ -64,6 +64,11 @@ function Base.show(io::IO, ::MIME"text/plain", A::LimitedMemoryMatrix) where {T}
     Base.showarg(io, A.value, true)
 end
 
+function Base.show(io::IO, A::LimitedMemoryMatrix) where {T}
+    print(io, Base.dims2string(size(A)), " ")
+    Base.showarg(io, A.value, true)
+end
+
 """
     hcat!(A::LimitedMemoryMatrix, B::AbstractVector)
 

test/lal.jl

@@ -238,35 +238,4 @@ end
 
         test_lal_identities(ld_results)
     end 
-end
-
-@testset "Cyclic Circulant Matrix" begin
-    rng = Random.Xoshiro(1234)
-    # 4-cyclic circulant matrix, creates blocks
-    I1 = Diagonal(fill(1.0, 3))
-    I2 = Diagonal(fill(1.0, 7))
-    I3 = Diagonal(fill(1.0, 3))
-    I4 = Diagonal(fill(1.0, 5))
-    B1 = rand(rng, size(I1, 1), size(I4, 2))
-    B2 = rand(rng, size(I2, 1), size(I1, 2))
-    B3 = rand(rng, size(I3, 1), size(I2, 2))
-    B4 = rand(rng, size(I4, 1), size(I3, 2))
-    Z12 = fill(0.0, size(I1, 1), size(I2, 2))
-    Z13 = fill(0.0, size(I1, 1), size(I3, 2))
-    Z14 = fill(0.0, size(I1, 1), size(I4, 2))
-    Z23 = fill(0.0, size(I2, 1), size(I3, 2))
-    Z24 = fill(0.0, size(I2, 1), size(I4, 2))
-    Z34 = fill(0.0, size(I3, 1), size(I4, 2))
-    Ac = [
-        I1   Z12  Z13 B1
-        B2   I2   Z23 Z24
-        Z13' B3   I3  Z34
-        Z14' Z24' B4  I4
-    ]
-    v = [rand(rng, size(I1, 1)); fill(0.0, size(Ac, 1) - size(I1, 1))]
-    w = [rand(rng, size(I1, 1)); fill(0.0, size(Ac, 1) - size(I1, 1))]
-    ld = IS.LookAheadLanczosDecomp(Ac, v, w; log=true, vw_normalized=false, max_block_size=10, max_memory=10, verbose=true)
-    ld_results = _iterate_and_collect_lal_intermediates(ld)
-    test_lal_identities(ld_results)
-end
-
+end