Add column pivoting to block_qr!

[leburgel]

I was trying to hack a way around #78 by using the BlockLanczos method to warm-start an eigsolve for 6 extremal eigenvectors using a previous partially converged result, and ran into what I think is a bug in block_qr! that actually made the warm start worse than starting from scratch.

block_qr! orthogonalized the vectors in the order they were given, so the detected rank could depend on that order. This adds column pivoting: the remaining vector with the largest norm is processed first.

This shows up when warm-starting BlockLanczos from the eigenvectors of an earlier, unconverged eigsolve. Those vectors share one residual direction, so AX - XΘ has rank 1 and the block should shrink to size 1 after the first step. The first vector is usually the best converged one, though. Normalizing its tiny residual first gives a direction dominated by rounding errors, and projecting the larger residuals onto it leaves remainders above qr_tol. The iteration then carries a spurious second vector, doubling the matvecs per step. Reversing the input order avoids it on master.

goodidx is now returned in pivot order, and R is upper triangular only up to a column permutation. The callers only rely on block[goodidx] * R reproducing the block, so they don't change.

Reproducer:

using KrylovKit, LinearAlgebra, Random
Random.seed!(1)
n = 4000
A = SymTridiagonal(2 .+ 0.5 .* rand(n), -ones(n - 1))
opts = (; tol = 1e-10, ishermitian = true, verbosity = 0, maxiter = 10_000)

# stop early, then warm-start BlockLanczos from the 6 returned eigenvectors
vals, vecs, info = eigsolve(A, rand(n), 6, :SR; opts..., krylovdim = 20, maxiter = 20)
println("normres = ", round.(info.normres; sigdigits = 2))
println("svdvals(AX - XΘ) = ", round.(svdvals(stack(A * v - λ * v for (λ, v) in zip(vals, vecs))); sigdigits = 2))
for X in (vecs, reverse(vecs)) # same vectors, reversed order
    iter = KrylovKit.BlockLanczosIterator(A, Block(copy.(X)), 100)
    fact = initialize(iter)
    expand!(iter, fact)
    println("block size after first expansion: ", fact.R_size)
end
for kd in (20, 60)
    cold = eigsolve(A, rand(n), 6, :SR; opts..., krylovdim = kd)[3].numops
    warm = eigsolve(A, Block(copy.(vecs)), 6, :SR; opts..., krylovdim = kd)[3].numops
    println("krylovdim = $kd: matvecs cold = $cold, warm = $warm")
end

master:

normres = [6.5e-9, 0.00047, 0.0041, 0.0096, 0.0073, 0.014]
svdvals(AX - XΘ) = [0.019, 5.0e-15, 4.5e-15, 2.6e-15, 2.1e-15, 1.7e-15]
block size after first expansion: 2
block size after first expansion: 1
krylovdim = 20: matvecs cold = 886, warm = 3589
krylovdim = 60: matvecs cold = 533, warm = 565

this PR:

normres = [6.5e-9, 0.00047, 0.0041, 0.0096, 0.0073, 0.014]
svdvals(AX - XΘ) = [0.019, 5.0e-15, 4.5e-15, 2.6e-15, 2.1e-15, 1.7e-15]
block size after first expansion: 1
block size after first expansion: 1
krylovdim = 20: matvecs cold = 886, warm = 1173
krylovdim = 60: matvecs cold = 533, warm = 397

Added a test in test/block.jl that fails on master.

With krylovdim = 20 the warm start is still slower than a cold start (1173 vs 886 matvecs), which was a bit disappointing since that was exactly the case I was actually interested in.Could this be a separate problem in the BlockLanczos restart?

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 86.69%. Comparing base (531b074) to head (5f0bbc2).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #170      +/-   ##
==========================================
+ Coverage   86.65%   86.69%   +0.04%     
==========================================
  Files          36       36              
  Lines        3918     3915       -3     
==========================================
- Hits         3395     3394       -1     
+ Misses        523      521       -2     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.