From 4c88a5908f71dc3f33697bceff53d5c8b8feaba7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olivier=20Gauth=C3=A9?= Date: Fri, 1 Aug 2025 16:22:03 -0400 Subject: [PATCH 1/6] working randomized SVD --- src/PEPSKit.jl | 1 + .../contractions/ctmrg_contractions.jl | 11 ++ src/algorithms/ctmrg/projectors.jl | 126 +++++++++++++++++- src/algorithms/ctmrg/simultaneous.jl | 27 +++- src/environments/ctmrg_environments.jl | 31 +++++ 5 files changed, 190 insertions(+), 6 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 9e1c24b74..bb41fcd13 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -5,6 +5,7 @@ using Compat using Accessors: @set, @reset using VectorInterface import VectorInterface as VI +import Random using TensorKit using TensorKit: TruncationScheme diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index 651d3a378..35df55f90 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -263,6 +263,17 @@ function contract_projectors(U, S, V, Q, Q_next) return P_left, P_right end +function contract_projectors(U, S, V, fq::FourQuadrants) + isqS = sdiag_pow(S, -0.5) + # use * to respect fermionic case + p1 = (codomainind(fq), domainind(fq)) + p2 = (codomainind(fq), (numout(fq) + 1,)) + P_left = tensorcontract(fq.Q3, p1, false, fq.Q4 * V' * isqS, p2, false, p2) + p3 = ((1,), codomainind(fq) .+ 1) + P_right = tensorcontract(isqS * U' * fq.Q1, p3, false, fq.Q2, p1, false, p3) + return P_left, P_right +end + """ half_infinite_environment(quadrant1, quadrant2) half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, A_1, A_2) diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl index f2dc247b3..fe77498f2 100644 --- a/src/algorithms/ctmrg/projectors.jl +++ b/src/algorithms/ctmrg/projectors.jl @@ -87,7 +87,7 @@ Construct the half-infinite projector algorithm based on the following keyword a - `:notrunc` : No singular values are truncated and the performed SVDs are exact - `:truncerr` : Additionally supply error threshold `η`; truncate to the maximal virtual dimension of `η` - `:truncdim` : Additionally supply truncation dimension `η`; truncate such that the 2-norm of the truncated values is smaller than `η` - - `:truncspace` : Additionally supply truncation space `η`; truncate according to the supplied vector space + - `:truncspace` : Additionally supply truncation space `η`; truncate according to the supplied vector space - `:truncbelow` : Additionally supply singular value cutoff `η`; truncate such that every retained singular value is larger than `η` * `verbosity::Int=$(Defaults.projector_verbosity)` : Projector output verbosity which can be: 0. Suppress output information @@ -125,7 +125,7 @@ Construct the full-infinite projector algorithm based on the following keyword a - `:notrunc` : No singular values are truncated and the performed SVDs are exact - `:truncerr` : Additionally supply error threshold `η`; truncate to the maximal virtual dimension of `η` - `:truncdim` : Additionally supply truncation dimension `η`; truncate such that the 2-norm of the truncated values is smaller than `η` - - `:truncspace` : Additionally supply truncation space `η`; truncate according to the supplied vector space + - `:truncspace` : Additionally supply truncation space `η`; truncate according to the supplied vector space - `:truncbelow` : Additionally supply singular value cutoff `η`; truncate such that every retained singular value is larger than `η` * `verbosity::Int=$(Defaults.projector_verbosity)` : Projector output verbosity which can be: 0. Suppress output information @@ -148,7 +148,7 @@ PROJECTOR_SYMBOLS[:fullinfinite] = FullInfiniteProjector Determine left and right projectors at the bond given determined by the enlarged corners and the given coordinate using the specified `alg`. """ -function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjector) +function compute_projector(enlarged_corners, coordinate, last_space, alg::HalfInfiniteProjector) # SVD half-infinite environment halfinf = half_infinite_environment(enlarged_corners...) svd_alg = svd_algorithm(alg, coordinate) @@ -166,7 +166,7 @@ function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjec P_left, P_right = contract_projectors(U, S, V, enlarged_corners...) return (P_left, P_right), (; U, S, V, info...) end -function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjector) +function compute_projector(enlarged_corners, coordinate, last_space, alg::FullInfiniteProjector) halfinf_left = half_infinite_environment(enlarged_corners[1], enlarged_corners[2]) halfinf_right = half_infinite_environment(enlarged_corners[3], enlarged_corners[4]) @@ -187,3 +187,121 @@ function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjec P_left, P_right = contract_projectors(U, S, V, halfinf_left, halfinf_right) return (P_left, P_right), (; U, S, V, info...) end + +# ========================================================================================== +# TBD: compute ncv vectors proj.ncv_ratio * old_nvectors? +struct RandomizedProjector{S, T, R} <: ProjectorAlgorithm + svd_alg::S + trscheme::T + rng::R + oversampling::Int64 + max_full::Int64 + verbosity::Int +end + +PROJECTOR_SYMBOLS[:randomized] = RandomizedProjector + +svd_algorithm(alg::RandomizedProjector) = alg.svd_alg + +#= +# TBD is this needed? +function RandomizedProjector(; + rng = Random.default_rng(), + svd_alg = (;), + trscheme = (;), + oversampling = 10, + max_full = 100, + verbosity = Defaults.projector_verbosity, + ) + + # parse SVD forward & rrule algorithm + svd_algorithm = _alg_or_nt(SVDAdjoint, svd_alg) + + # parse truncation scheme + truncation_scheme = if trscheme isa TruncationScheme + trscheme + elseif trscheme isa NamedTuple + _TruncationScheme(; trscheme...) + else + throw(ArgumentError("unknown trscheme $trscheme")) + end + + return RandomizedProjector( + svd_algorithm, truncation_scheme, rng, oversampling, max_full, verbosity + ) +end +=# + + +# needed as default interface in PEPSKit.ProjectorAlgorithm +function RandomizedProjector(svd_algorithm, trscheme, verbosity) + @show "Hi RandomizedProjector" + @show which( + RandomizedProjector, typeof.( + ( + svd_algorithm, trscheme, Random.default_rng(), 10, 100, verbosity, + ) + ) + ) + return RandomizedProjector( + svd_algorithm, trscheme, Random.default_rng(), 10, 100, verbosity + ) +end + +function random_domain(alg::RandomizedProjector, full_space, last_space) + sector_dims = map(sectors(full_space)) do s + if dim(full_space, s) <= alg.max_full + n = dim(full_space, s) + else + n = dim(last_space, s) + alg.oversampling + end + return s => n + end + return Vect[sectortype(last_space)](sector_dims) +end + + +function randomized_range_finder(A::AbstractTensorMap, alg::RandomizedProjector, randomized_space) + Q = TensorMap{eltype(A)}(undef, domain(A) ← randomized_space) + foreach(blocks(Q)) do (s, b) + m, n = size(b) + if m <= alg.max_full + b .= LinearAlgebra.I(m) + else + Ω = randn(alg.rng, eltype(A), domain(A) ← Vect[sectortype(A)](s => n)) + Y = A * Ω + Qs, _ = leftorth!(Y) + b .= block(Qs, s) + end + end + return Q +end + +# impose full env, could also be defined for half_infinite_environment, little gain +function compute_projector(fq, coordinate, last_space, alg::RandomizedProjector) + full_space = fuse(domain(fq)) + randomized_space = random_domain(alg, full_space, last_space) + Q = randomized_range_finder(fq, alg, randomized_space) + B = Q' * fq + + svd_alg = svd_algorithm(alg, coordinate) + U′, S, V, info = tsvd!(B, svd_alg; trunc = alg.trscheme) + U = Q * U′ + foreach(blocks(S)) do (s, b) + if size(b, 1) == dim(randomized_space, s) && size(b, 1) < dim(full_space, s) + @warn("Sector is too small, kept all computed values: ", s) + end + end + + # Check for degenerate singular values; still needed for exact blocks + Zygote.isderiving() && ignore_derivatives() do + if alg.verbosity > 0 && is_degenerate_spectrum(S) + svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S)) + @warn("degenerate singular values detected: ", svals) + end + end + + @reset info.truncation_error = info.truncation_error / norm(S) # normalize truncation error + P_left, P_right = contract_projectors(U, S, V, fq) + return (P_left, P_right), (; U, S, V, info...) +end diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index ffe68bec2..ec0eb12c8 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -97,7 +97,8 @@ function simultaneous_projectors( trscheme = truncation_scheme(alg, env.edges[coordinate[1], coordinate′[2:3]...]) alg′ = @set alg.trscheme = trscheme ec = (enlarged_corners[coordinate...], enlarged_corners[coordinate′...]) - return compute_projector(ec, coordinate, alg′) + last_space = space(env.edges[coordinate[1], coordinate′[2:3]...], 1) + return compute_projector(ec, coordinate, last_space, alg′) end function simultaneous_projectors( coordinate, enlarged_corners::Array{E, 3}, env, alg::FullInfiniteProjector @@ -109,13 +110,35 @@ function simultaneous_projectors( coordinate2 = _next_coordinate(coordinate, rowsize, colsize) coordinate3 = _next_coordinate(coordinate2, rowsize, colsize) coordinate4 = _next_coordinate(coordinate3, rowsize, colsize) + last_space = space(env.edges[coordinate[1], coordinate′[2:3]...], 1) ec = ( enlarged_corners[coordinate4...], enlarged_corners[coordinate...], enlarged_corners[coordinate2...], enlarged_corners[coordinate3...], ) - return compute_projector(ec, coordinate, alg′) + return compute_projector(ec, coordinate, last_space, alg′) +end + +# TBD share code with FullInfiniteProjector? +function simultaneous_projectors( + coordinate, enlarged_corners::Array{E, 3}, env, alg::RandomizedProjector + ) where {E} + coordinate′ = _next_coordinate(coordinate, size(env)[2:3]...) + trscheme = truncation_scheme(alg, env.edges[coordinate[1], coordinate′[2:3]...]) + alg′ = @set alg.trscheme = trscheme + rowsize, colsize = size(enlarged_corners)[2:3] + coordinate2 = _next_coordinate(coordinate, rowsize, colsize) + coordinate3 = _next_coordinate(coordinate2, rowsize, colsize) + coordinate4 = _next_coordinate(coordinate3, rowsize, colsize) + fq = FourQuadrants( + enlarged_corners[coordinate4...], + enlarged_corners[coordinate...], + enlarged_corners[coordinate2...], + enlarged_corners[coordinate3...], + ) + last_space = space(env.edges[coordinate[1], coordinate′[2:3]...], 1) + return compute_projector(fq, coordinate, last_space, alg′) end """ diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index ef269d270..4a2e9b9f7 100644 --- a/src/environments/ctmrg_environments.jl +++ b/src/environments/ctmrg_environments.jl @@ -477,3 +477,34 @@ function VI.inner(env₁::CTMRGEnv, env₂::CTMRGEnv) return inner((env₁.corners, env₁.edges), (env₂.corners, env₂.edges)) end VI.norm(env::CTMRGEnv) = norm((env.corners, env.edges)) + +# ========================================================================================== + +struct FourQuadrants{S, T, N, TM} <: AbstractTensorMap{S, T, N, N} + Q1::TM + Q2::TM + Q3::TM + Q4::TM + + function FourQuadrants{T, S, N, TM}( + Q1, Q2, Q3, Q4 + ) where {S, T, N, TM <: AbstractTensorMap{T, S, N, N}} + return new{T, S, N, TM}(Q1, Q2, Q3, Q4) + end +end + +function FourQuadrants(Q1, Q2, Q3, Q4) + return FourQuadrants{eltype(Q1), spacetype(Q1), numin(Q1), typeof(Q1)}(Q1, Q2, Q3, Q4) +end + +TensorKit.TensorMap(fq::FourQuadrants) = fq.Q1 * fq.Q2 * fq.Q3 * fq.Q4 + +TensorKit.space(fq::FourQuadrants) = codomain(fq.Q1) ← domain(fq.Q4) + +# TBD use tensorcontract to handle fermions? +(fq::FourQuadrants)(m) = fq.Q1 * (fq.Q2 * (fq.Q3 * (fq.Q4 * m))) + +Base.:*(m::AbstractTensorMap, fq::FourQuadrants) = m * fq.Q1 * fq.Q2 * fq.Q3 * fq.Q4 +Base.:*(fq::FourQuadrants, m::AbstractTensorMap) = fq.Q1 * (fq.Q2 * (fq.Q3 * (fq.Q4 * m))) + +Base.adjoint(fq::FourQuadrants) = FourQuadrants(fq.Q4', fq.Q3', fq.Q2', fq.Q1') From 05d078ed2aa004844ab8f70df03ea0bc0a1d1e04 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olivier=20Gauth=C3=A9?= Date: Thu, 18 Sep 2025 19:57:49 -0400 Subject: [PATCH 2/6] cleaning --- src/algorithms/ctmrg/projectors.jl | 41 +++--------------------------- 1 file changed, 3 insertions(+), 38 deletions(-) diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl index fe77498f2..727e32a9c 100644 --- a/src/algorithms/ctmrg/projectors.jl +++ b/src/algorithms/ctmrg/projectors.jl @@ -203,48 +203,13 @@ PROJECTOR_SYMBOLS[:randomized] = RandomizedProjector svd_algorithm(alg::RandomizedProjector) = alg.svd_alg -#= -# TBD is this needed? -function RandomizedProjector(; - rng = Random.default_rng(), - svd_alg = (;), - trscheme = (;), - oversampling = 10, - max_full = 100, - verbosity = Defaults.projector_verbosity, - ) - - # parse SVD forward & rrule algorithm - svd_algorithm = _alg_or_nt(SVDAdjoint, svd_alg) - - # parse truncation scheme - truncation_scheme = if trscheme isa TruncationScheme - trscheme - elseif trscheme isa NamedTuple - _TruncationScheme(; trscheme...) - else - throw(ArgumentError("unknown trscheme $trscheme")) - end - - return RandomizedProjector( - svd_algorithm, truncation_scheme, rng, oversampling, max_full, verbosity - ) -end -=# - +_default_randomized_oversampling = 10 +_default_randomized_max_full = 100 # needed as default interface in PEPSKit.ProjectorAlgorithm function RandomizedProjector(svd_algorithm, trscheme, verbosity) - @show "Hi RandomizedProjector" - @show which( - RandomizedProjector, typeof.( - ( - svd_algorithm, trscheme, Random.default_rng(), 10, 100, verbosity, - ) - ) - ) return RandomizedProjector( - svd_algorithm, trscheme, Random.default_rng(), 10, 100, verbosity + svd_algorithm, trscheme, Random.default_rng(), _default_randomized_oversampling, _default_randomized_max_full, verbosity ) end From cd1a1c3a773aab2eaeca43a796950f0c2f2adfa1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olivier=20Gauth=C3=A9?= Date: Fri, 19 Sep 2025 18:13:36 -0400 Subject: [PATCH 3/6] add subspace iterations --- src/algorithms/ctmrg/projectors.jl | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl index 727e32a9c..853b91292 100644 --- a/src/algorithms/ctmrg/projectors.jl +++ b/src/algorithms/ctmrg/projectors.jl @@ -194,8 +194,9 @@ struct RandomizedProjector{S, T, R} <: ProjectorAlgorithm svd_alg::S trscheme::T rng::R - oversampling::Int64 - max_full::Int64 + oversampling::Int + max_full::Int + n_subspace_iter::Int verbosity::Int end @@ -205,11 +206,13 @@ svd_algorithm(alg::RandomizedProjector) = alg.svd_alg _default_randomized_oversampling = 10 _default_randomized_max_full = 100 +_default_n_subspace_iter = 2 # needed as default interface in PEPSKit.ProjectorAlgorithm function RandomizedProjector(svd_algorithm, trscheme, verbosity) return RandomizedProjector( - svd_algorithm, trscheme, Random.default_rng(), _default_randomized_oversampling, _default_randomized_max_full, verbosity + svd_algorithm, trscheme, Random.default_rng(), _default_randomized_oversampling, + _default_randomized_max_full, _default_n_subspace_iter, verbosity ) end @@ -229,11 +232,17 @@ end function randomized_range_finder(A::AbstractTensorMap, alg::RandomizedProjector, randomized_space) Q = TensorMap{eltype(A)}(undef, domain(A) ← randomized_space) foreach(blocks(Q)) do (s, b) + Aad = A' m, n = size(b) if m <= alg.max_full b .= LinearAlgebra.I(m) else Ω = randn(alg.rng, eltype(A), domain(A) ← Vect[sectortype(A)](s => n)) + for _ in 1:alg.n_subspace_iter + Y = A * Ω + Qs, _ = leftorth!(Y) + Ω, _ = leftorth!(Aad * Qs) + end Y = A * Ω Qs, _ = leftorth!(Y) b .= block(Qs, s) @@ -248,6 +257,7 @@ function compute_projector(fq, coordinate, last_space, alg::RandomizedProjector) randomized_space = random_domain(alg, full_space, last_space) Q = randomized_range_finder(fq, alg, randomized_space) B = Q' * fq + normalize!(B) # TODO better way? svd_alg = svd_algorithm(alg, coordinate) U′, S, V, info = tsvd!(B, svd_alg; trunc = alg.trscheme) From 3b0c96af8c4cb3197daaf6dd44265c39ab12708a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olivier=20Gauth=C3=A9?= Date: Wed, 15 Oct 2025 14:21:10 +0200 Subject: [PATCH 4/6] quickfix for SequentialProjector --- src/algorithms/ctmrg/sequential.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index 04be3285b..ecf742bd4 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -101,7 +101,7 @@ function sequential_projectors( r′ = _prev(r, size(env, 2)) Q1 = TensorMap(EnlargedCorner(network, env, (SOUTHWEST, r, c))) Q2 = TensorMap(EnlargedCorner(network, env, (NORTHWEST, r′, c))) - return compute_projector((Q1, Q2), coordinate, alg) + return compute_projector((Q1, Q2), coordinate, nothing, alg) end function sequential_projectors( coordinate::NTuple{3, Int}, network, env::CTMRGEnv, alg::FullInfiniteProjector @@ -116,7 +116,7 @@ function sequential_projectors( TensorMap(EnlargedCorner(network, env, coordinate_nw)), TensorMap(EnlargedCorner(network, env, coordinate_ne)), ) - return compute_projector(ec, coordinate, alg) + return compute_projector(ec, coordinate, nothing, alg) end """ From a2a99437061ebdc9fbfb9486b2c440afa35e6fa5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olivier=20Gauth=C3=A9?= Date: Thu, 13 Nov 2025 13:34:32 +0100 Subject: [PATCH 5/6] update to MatrixAlgebraKit --- src/algorithms/ctmrg/projectors.jl | 8 ++++---- src/algorithms/ctmrg/simultaneous.jl | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl index 95893cfa9..5b5672fe2 100644 --- a/src/algorithms/ctmrg/projectors.jl +++ b/src/algorithms/ctmrg/projectors.jl @@ -193,7 +193,7 @@ end # TBD: compute ncv vectors proj.ncv_ratio * old_nvectors? struct RandomizedProjector{S, T, R} <: ProjectorAlgorithm svd_alg::S - trscheme::T + trunc::T rng::R oversampling::Int max_full::Int @@ -210,9 +210,9 @@ _default_randomized_max_full = 100 _default_n_subspace_iter = 2 # needed as default interface in PEPSKit.ProjectorAlgorithm -function RandomizedProjector(svd_algorithm, trscheme, verbosity) +function RandomizedProjector(svd_algorithm, trunc, verbosity) return RandomizedProjector( - svd_algorithm, trscheme, Random.default_rng(), _default_randomized_oversampling, + svd_algorithm, trunc, Random.default_rng(), _default_randomized_oversampling, _default_randomized_max_full, _default_n_subspace_iter, verbosity ) end @@ -261,7 +261,7 @@ function compute_projector(fq, coordinate, last_space, alg::RandomizedProjector) normalize!(B) # TODO better way? svd_alg = svd_algorithm(alg, coordinate) - U′, S, V, info = tsvd!(B, svd_alg; trunc = alg.trscheme) + U′, S, V, info = svd_trunc!(B, svd_alg; trunc = alg.trunc) U = Q * U′ foreach(blocks(S)) do (s, b) if size(b, 1) == dim(randomized_space, s) && size(b, 1) < dim(full_space, s) diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index 68a8a1e45..af8a9cbe5 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -125,8 +125,8 @@ function simultaneous_projectors( coordinate, enlarged_corners::Array{E, 3}, env, alg::RandomizedProjector ) where {E} coordinate′ = _next_coordinate(coordinate, size(env)[2:3]...) - trscheme = truncation_scheme(alg, env.edges[coordinate[1], coordinate′[2:3]...]) - alg′ = @set alg.trscheme = trscheme + trunc = truncation_strategy(alg, env.edges[coordinate[1], coordinate′[2:3]...]) + alg′ = @set alg.trunc = trunc rowsize, colsize = size(enlarged_corners)[2:3] coordinate2 = _next_coordinate(coordinate, rowsize, colsize) coordinate3 = _next_coordinate(coordinate2, rowsize, colsize) From 27b9a683412dfbabae79fae4a8399c0aaa4fd942 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olivier=20Gauth=C3=A9?= Date: Mon, 5 Jan 2026 12:02:20 +0100 Subject: [PATCH 6/6] adapt to MatrixAlgebraKit left_orth! --- src/algorithms/ctmrg/projectors.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl index 5b5672fe2..ac346acfb 100644 --- a/src/algorithms/ctmrg/projectors.jl +++ b/src/algorithms/ctmrg/projectors.jl @@ -241,11 +241,11 @@ function randomized_range_finder(A::AbstractTensorMap, alg::RandomizedProjector, Ω = randn(alg.rng, eltype(A), domain(A) ← Vect[sectortype(A)](s => n)) for _ in 1:alg.n_subspace_iter Y = A * Ω - Qs, _ = leftorth!(Y) - Ω, _ = leftorth!(Aad * Qs) + Qs, _ = left_orth!(Y) + Ω, _ = left_orth!(Aad * Qs) end Y = A * Ω - Qs, _ = leftorth!(Y) + Qs, _ = left_orth!(Y) b .= block(Qs, s) end end