Nonlinear and linear solver options
Methods, preconditioners, refresh policies, and factorization choices.
See the API overview for other analyses and conventions for units and array ordering.
JosephsonCircuits.Always Type
Always()Rebuild the preconditioner before every Newton step: the default, and the policy whose solve path depends on nothing measured, so that the same problem solves the same way every time. Probe measures whether a rebuild pays and skips the ones which do not.
JosephsonCircuits.Automatic Type
Automatic()The preconditioner chosen when it is built, by the number of tones and the memory the factors would take:
one tone:
FullJacobianwith the backend's sparse factorization, whatever the memory;two or more tones:
FullJacobianwith a single precisionBlockFactorizationwhen its factors, sized exactly from the symbolic analysis byblockfactorbytes, take at most half the free memory of the backend (freememory), leaving the rest to the system, the Krylov basis and the products; otherwiseMeasuredBandwith the backend's sparse factorization.
The choice looks at nothing else, neither the mixing order of the tones nor the circuit's topology, so that it generalizes. On one tone the full Jacobian's sparse factors give one Krylov iteration per Newton step where the block diagonal needs many, and no escalation; they are taken without sizing them first, since they grow with the square of the mode count and one tone keeps it low, and a problem too large for them wants an explicit BlockDiagonal. See resolveautomatic.
JosephsonCircuits.Backtracking Type
Backtracking(; interpolate = true, safeguardlow = 0.1,
safeguardhigh = 0.5, c1 = 1e-4, maxbacktracks = 10, maxfailures = 2)Armijo backtracking along the step, the line search of every method. A trial which meets the Armijo condition with constant c1 is accepted; otherwise the step is shortened and retried, at most maxbacktracks times, and maxfailures consecutive steps which return the best decreasing trial rather than an Armijo step are a stall.
interpolate chooses how the step is shortened. Interpolating fits the merit function along the step, a quadratic through the full step for the first proposal and a cubic through the two latest trials for each one after, and clamps every proposal to [safeguardlow, safeguardhigh] of the trial before it; it suits an exact step, whose usable length is a particular number. Halving makes no fit: the first backtrack is to half the step and each later one multiplies it by safeguardhigh, so only that safeguard acts; it suits an inexact step, which is perturbed enough that a power of one half lands acceptably. Halving ends a hopeless solve by itself, since its coarse trials soon fail the Armijo condition, where interpolation always finds a step short enough to meet it and leaves the ending to the stall rule (residualstalled).
safeguardlow bounds how far one interpolated trial may cut. A Newton step in NewtonKrylov refreshes the preconditioner and a backtrack does not, so too deep a cut buys fewer trials at the price of more steps.
JosephsonCircuits.BlockFactorization Type
BlockFactorization(singletons = nothing; precision = nothing, refine = true)The AbstractFactorization of a ModeCouplingPreconditioner by dense blocks over the circuit graph rather than by a scalar sparse solver, and, handed a sparse matrix with a block size, of the linearized solve.
The harmonic balance Jacobian has two structures: its sparsity follows the circuit graph, and every nonlinear connection carries a dense coupling between harmonic modes. A scalar sparse solver sees only their product and rediscovers the circuit's block structure, entry by entry, in its symbolic analysis. This factorization eliminates the circuit graph instead, treating every circuit node as a supernode whose block holds the real-layout slots of all retained modes at that node: the nodes are ordered by KLU's analysis of the circuit-node graph (block triangular form, then a fill reducing ordering), merged along the elimination tree into larger supernodes only where the merged block holds no entry the factors would not hold node by node (the separators of a meshed circuit; a chain stays node by node, see amalgamate), and eliminated with pivoted dense LU on the diagonal blocks, dense products for the Schur updates and dense matrix-vector products for the solves. The blocks are assembled on the backend straight from the Fourier coefficients with realstructureentry, the same per entry value the sparse assembly uses, so no sparse Jacobian is ever formed. The structure comes from the circuit graph alone, whatever the circuit is.
The same specification serves two operators. Handed a BlockJacobian it factorizes the preconditioner's mode clusters over the circuit graph, as above; handed a SparseMatrixCSC with a blocksize (see factorize) it builds a SparseBlockFactorization, the direct solve of the linearized system in dense node blocks, batched over the frequencies of a device sweep, with precision then meaning the precision of the factors of a double matrix (equilibrated and refined when single). That solve merges no nodes: a merged supernode of a meshed circuit can be ill conditioned at a frequency where the matrix is not, and pivoting within it then costs the solution digits (blocksymbolic), where in the preconditioner it costs only iterations.
The linearized solve takes it at two or more tones, on the host and on a device, when its factors fit the memory budget (linearizedfactorization): several tones make the node blocks large and dense, and BLAS-3 kernels eliminate them, and the separators a junction lattice fills, faster than a scalar factorization does. A scalar factorization, KLUfactorization on the host, pivots across the whole matrix where this one pivots within a supernode, and keeps the sparsity within the node blocks which this one stores dense.
The coupling set of the preconditioner is honored at the level of its clusters: the retained coupling graph of the modes is split into its connected components, every component of two or more modes becomes one block factorization over the circuit graph restricted to those modes' slots, and the modes left single are solved by the mode block diagonal. FullJacobian is therefore one factorization of the complete Jacobian, an exact solve; a mask made of complete clusters (as Clusters produces) is one factorization per cluster; and a coupling set which is not a union of complete clusters, a HarmonicBand say, is factorized on its closure, which keeps at least every coupling the set asked for.
Factor storage grows as the square of the number of retained slots per node and the arithmetic as its cube, which is what bounds it.
singletons is the sparse AbstractFactorization of the block diagonal of the modes left single, the backend's default (KLU on the host, cuDSS on a device) when nothing. precision is the floating point type of the blocks, that of the preconditioner when nothing: Float32 is the mixed precision form, which halves the storage and runs at a device's single precision rate while the iteration stays in double precision.
The block factorization does not pivot across its supernodes. A supernode whose diagonal block is singular, which happens when the part of the circuit eliminated up to it resonates at a mode frequency with the rest held at zero, stops it with a SingularException, though the matrix is not singular. The mode coupling preconditioner then falls back to the backend's sparse factorization (see refactorize!); the linearized solve does too when it chose the block factorization itself, solving its sweep again with KLUfactorization() or CUDSSFactorization(), and throws when it was given one.
refine concerns the linearized solve: single precision factors of a double system refine their solutions against the double residual to double accuracy by default; refine = false leaves them single precision solutions computed entirely in single precision, equilibrated, with no refinement, for the cases where single precision scattering parameters suffice, at the speed of a device whose single precision rate far exceeds its double one; the accuracy then falls with the conditioning of the system, furthest on strongly resonant multi-tone lines. The outputs are returned in double either way.
JosephsonCircuits.CUDSSFactorization Type
CUDSSFactorization(; precision = nothing, kwargs...)An AbstractFactorization backed by NVIDIA's cuDSS direct sparse solver, for use on a GPU.
precision is the floating point type the factors are held and the triangular solves run in, nothing for the iteration's own. It is a setting of a NewtonKrylov preconditioner: the factors of a preconditioner need only make the Krylov solve converge, so Float32 halves their memory and runs them at a device's single precision rate, the preconditioner assembling its matrix in that precision and converting the residual and the correction around each solve. Such a preconditioner is not exact, and a Krylov solve it fails is escalated to the same factors in the iteration's precision (escalatepreconditioner!). A direct solve factorizes in the precision of its iteration and refuses a factorization asking for another. cuDSS and BlockFactorization honour a precision; KLU and UMFPACK factorize in double whatever they are handed.
Requires CUDSS.jl and CUDA.jl to be loaded; without them the returned factorization raises an informative error when used, so the constructor itself is always available and a script can select it unconditionally. Nothing else in the solver changes: AbstractFactorization already separates factorize (the symbolic analysis and the first numeric factorization) from refactorize! (a numeric refactorization reusing the analysis), which is exactly the split cuDSS wants.
That split is what makes this worth doing. The symbolic analysis depends only on the circuit topology and the retained mode coupling, neither of which changes across Newton steps, while the values change at every step; the analysis is almost all of the cost of a factorization and is paid once.
cuDSS is handed the whole matrix, block diagonal or not: it discovers the independent blocks of a mode block diagonal itself and works on them together, which is faster than handing it the blocks as a uniform batch.
sourceJosephsonCircuits.ExternalSolver Type
ExternalSolver(f)Solve the operating point with a caller supplied root finder.
f(prob, u0) receives an HBNonlinearProblem and the initial value in the real representation, and returns (u, converged). Everything it needs is on prob: hbresidual!, hbjvp!, JacobianOperator, preconditioner, and the tolerance of the solve, prob.atol. The root is held to that tolerance as well: a u whose residual norm exceeds it has not converged, whatever converged says.
This is the plug point for a solver the package does not know about. A NonlinearSolve.jl algorithm, a hand written continuation stepper or a homotopy all go here without an extension.
The assembled real Jacobian is available on the problem unless assemblejacobian = false was passed, which is what a matrix-free solver wants: on a multi-tone problem that plan is the largest object in the solve. The problem lives on the host, so a solve with an external solver runs on the default backend = CPU().
ExternalSolver() do prob, u0
u = copy(u0); F = similar(u)
hbresidual!(F, prob, u)
# built once, and refactorized at each new point on its structure
P = preconditioner(prob, u)
for k in 1:40
norm(F) <= prob.atol && return (u, true)
J = JacobianOperator(prob, u)
d, st = Krylov.gmres(J, -F; N = P, rtol = 1e-10, atol = 0.0)
st.solved || return (u, false)
u .+= d; hbresidual!(F, prob, u)
JosephsonCircuits.updatepreconditioner!(P, u)
end
return (u, norm(F) <= prob.atol)
endAbsolute tolerances stall Newton
Note atol = 0.0. Krylov.jl defaults to atol = sqrt(eps()), about 1.5e-8, and stops as soon as the linear residual falls below it. Once the Newton residual is smaller than that – which is the whole point of the last few Newton steps – every linear solve returns immediately having done zero iterations, reports success, and hands back a zero step. Newton then stagnates while nothing reports a failure. Any external Krylov solver used inside a Newton loop wants its absolute tolerance set to zero and its stopping left to the relative one.
JosephsonCircuits.GMRES Type
GMRES(; restart = 400, maxrestarts = 4)The restarted GMRES of this package, with Givens rotations, the recycling subspace harvest and the preconditioner escalation the solver was built around. The default, and the only solver supporting deflation recycling, because the harvest (harvestcycle!) reads the Arnoldi basis out of the internal workspace. restart is the cycle length and maxrestarts the restart budget per solve. A long cycle is the default because a restricted preconditioner leaves a few directions a short Krylov space cannot resolve, and a restart discards the progress on them; the basis of restart + 1 vectors is cheap next to the sparse factorization an escalation would build.
JosephsonCircuits.KLUfactorization Type
KLUfactorization(; pivottol = 1e-6, kwargs...)The AbstractFactorization using KLU.jl, a sparse LU factorization suited to circuit matrices. This is the default on the host. The fill reducing ordering is chosen by kluordered rather than left at KLU's default.
A refactorization reuses the pivot sequence of the factorization before it, which the values of a later matrix can make unstable: a pivot which has become small inflates the entries its elimination updates, and the solution loses about as many digits as they grow. pivottol bounds that growth. When the largest pivot of the refactorized matrix, whose rows KLU scales to a largest entry of one, exceeds 1/pivottol, the same values are factorized again with fresh partial pivoting, reusing the symbolic analysis (klurefactor!). Partial pivoting keeps the growth of a fresh factorization far below the default bound, so a stable refactorization is not repeated; pivottol = 0 keeps every refactorization. kwargs are passed to KLU.klu.
JosephsonCircuits.KrylovJL Type
KrylovJL(method::Symbol = :gmres; kwargs...)Solve the Newton step with Krylov.jl: :gmres, :fgmres, :bicgstab, :dqgmres or any other solver taking an operator and a right hand side. Requires Krylov.jl to be loaded; the method lives in the package extension.
Only the linear solve changes. The forcing term, the line search, the preconditioner escalation and the stagnation handling are untouched, and the mode coupling preconditioner is passed through unchanged because it is applied by mul!, which is what Krylov.jl's N argument consumes. Deflation recycling is unavailable, since it depends on the internal workspace. The Krylov.jl workspace, which allocates a restarting method's whole basis when it is made, is made at the first linear solve of a system and kept for the rest; kwargs naming its size (memory, window) go to it, the others to each solve.
JosephsonCircuits.LUfactorization Type
LUfactorization(; kwargs...)The AbstractFactorization using the UMFPACK sparse LU factorization LinearAlgebra.lu, with kwargs passed to it.
JosephsonCircuits.MeasuredBand Type
MeasuredBand(; tol = 1e-2, budget = 0.25, factorization = nothing)A HarmonicBand whose per tone width is measured from the Fourier coefficients of cos(phi(t)) at every point and widened when the drive demands it, starting from the block diagonal; see cosphibandwidths for tol and budget. For two strong tones this is the setting that matters.
JosephsonCircuits.Never Type
Never()Rebuild the preconditioner only when it is forced: a linear solve which makes progress but misses its tolerance, a direction which is not a descent direction, a line search which finds no decrease where the rebuild can change the step, a successful escalation, or the recovery of a residual history which has stalled (residualstalled), which takes exact Newton steps from a rebuilt preconditioner. A stagnated solve is not retried and does not rebuild; its step is replaced by the preconditioner solve. The slow-solve report to the preconditioner (stalled!) is off as well, so a Clusters preconditioner never remeasures under this policy. A frozen preconditioner, for when a deflation (Floquet) is to carry the solve across the Newton path against a base built once.
JosephsonCircuits.Newton Type
Newton(; factorization = nothing, linesearch = Backtracking())Newton's method on the equivalent real system with the exact assembled real Jacobian, factorized by factorization (when nothing, KLU on the host and cuDSS on a device), the length of every step chosen by linesearch, a Backtracking, interpolating by default.
JosephsonCircuits.NewtonKrylov Type
NewtonKrylov(; preconditioner = Automatic(), linearsolver = GMRES(),
refresh = Always(), escalate = true, linesearch = Backtracking(),
precision = Float64)Jacobian-free Newton-Krylov with the mode coupling preconditioner: the default. preconditioner is an AbstractPreconditionerSpec; the default Automatic picks the full Jacobian for one tone and, for more, the set whose factors fit in memory. linearsolver is a GMRES or a KrylovJL solver, refresh Always (the default), Probe (which rebuilds the preconditioner only when a measured probe says a rebuild pays, at the price of a solve path which depends on measured times and so can differ between two runs) or Never. escalate allows a preconditioner which fails to reach its tolerance to be grown, within the memory the grown factors are predicted to take: the full set with its factors in less precision than the iteration, as the default single precision block factors of two or more tones are, to the iteration's precision; a band by one offset per tone; any other set to the full Jacobian (see escalatepreconditioner!). A refused escalation is recorded and the solve carries on.
The solve ends promptly when it cannot succeed and says why, in the reason of its IterationInfo: :iterations when the Newton steps are spent; :work when the Arnoldi steps exceed iterations restart lengths, so that a preconditioner which runs every linear solve to its limit cannot turn the step budget into hours; :linesearch when no sufficient decrease can be found (a step with no decrease at all is retried once from a rebuilt preconditioner when the rebuild can change it, see nlsolvekrylov!, and ends the solve otherwise or if it fails again, as does a direction which is still not a descent direction after the exact rescue, or two consecutive steps short of the Armijo condition); :progress when the residual history stopped coming down or comes down too slowly to reach the tolerance within the remaining budget, its rate is not improving and it is not accelerating, after one recovery which rebuilds the preconditioner and takes exact Newton steps from then on (residualstalled); :nonfinite when the residual norm at the initial point is not finite. A stall outside the Newton basin is the continuation problem Staged exists for. linesearch is the Backtracking which chooses the length of every step, interpolating by default; halving (Backtracking(interpolate = false)) suits an inexact preconditioner such as BlockDiagonal. precision is the floating point type of the iteration: the system on the backend, the Krylov vectors, and the matrix the preconditioner assembles and factorizes, except that KLU and UMFPACK factorize in double precision whatever they are handed, and a CUDSSFactorization or a BlockFactorization with a precision of its own holds its factors in that (factorizationprecision). The residual tolerance is raised to the rounding floor of the source in that precision, as atol describes.
The forcing sequence (Eisenstat-Walker choice 2 clamped to [1e-10, 0.9], starting at 0.3), the thresholds of a stagnated and of a slow linear solve and the floor of the linear solves are keywords of nlsolvekrylov!, with these defaults, rather than options of the method.
JosephsonCircuits.Probe Type
Probe()Decide each rebuild by measurement: one application of the stale preconditioner to the residual and one product give the one-step reduction rho = |J P^-1 F - F|/|F|; the same measurement on the fresh preconditioner, with the Arnoldi count k_fresh of its solve, calibrates the prediction k = k_fresh log(rho_fresh)/log(rho) of the stale solve's Arnoldi count, and the rebuild is skipped when k steps at the measured cost of a step are cheaper than the measured rebuild plus a fresh solve. Everything is measured, so the rule adapts to the device and the factorization; it pays when a rebuild is expensive next to a solve, as with a BlockFactorization of three tones. A rebuild forced by a failed, stalled or non-descent solve is never skipped. An exact preconditioner, the full coupling set factorized in the iteration's precision, is rebuilt at every step as under Always and not probed: its fresh reduction is roundoff, against which any stale one predicts a rebuild. Because the decision rests on measured times, the path a solve takes, and the answer within the tolerance, can differ between two runs of the same problem; the default Always is reproducible.
JosephsonCircuits.QRfactorization Type
QRfactorization(; kwargs...)The AbstractFactorization using the SPQR sparse QR factorization LinearAlgebra.qr, with kwargs passed to it. QR does not support refactorization in place, so each call factorizes from scratch.
JosephsonCircuits.QuasiNewton Type
QuasiNewton(; anderson = 5, factorization = nothing,
linesearch = Backtracking())The holomorphic Jacobian approximation with Anderson acceleration of depth anderson (the maximum number of previous iterates used for the extrapolation; zero disables it), factorized by factorization (when nothing, KLU on the host and cuDSS on a device), the length of every step chosen by linesearch, a Backtracking, interpolating by default; the curvilinear search of a rejected Anderson candidate follows it too. The harmonic balance residual is not complex differentiable, so this Jacobian is an approximation.
The zero frequency flux is complex
This method solves for a complex flux at every mode, the zero frequency mode included, so the imaginary part of a node's zero frequency flux is not held at zero and the converged value can carry a spurious imaginary part. The real part is the direct current flux; read only that, and prefer Newton or NewtonKrylov, which solve the real system, for a circuit with a direct current bias.
JosephsonCircuits.Staged Type
Staged(; grids = nothing, s0 = 0.5, smin = 0.02, interioratol = 1e-7,
interioriterations = 60, inner = NewtonKrylov(),
interiorescalation = false, maxattempts = 60, verbose = false)Source continuation on an adaptively grown harmonic grid, with inner solving every stage; see stagedhbnlsolve for the keywords. grids is the ladder of retained harmonic caps, nothing for the default ladder of the problem's Nharmonics.