docs menu

Module Jera.Linear

Linear systems.

A system is a linear function a on values of a structure 'x and a right-hand side r: its solution is the u with a u = r. The float tensors of a value are its vector, of one dtype, the solve's; its other tensors are carried from r. A solver says how a is used: dense materialises it from its products and factors it, banded does so for a matrix with entries near its diagonal only, and cg and gmres iterate on its products.

ocaml

(* u with (k I + L) u = r, L a matrix given by its product *)
let s =
  Linear.solve Nx.Ptree.tensor Linear.dense
    (fun u -> Nx.add (Nx.mul_s u k) (Nx.matmul l u))
    r

Check. Every answer is checked by one more application of a: its residual a u − r must meet the solver's bound, stated with each. A miss, or a non-finite u from a finite system, ends the lane Stalled: a is not linear, or it is singular or too ill-conditioned for the solver. A non-finite r or product of a ends it Not_finite. The answer's error is the magnitude of its residual, per component.

Derivative. The answer is stated as the zero of a u − r through Rune.root, with the solver as its linear solve: a derivative solves a system of a, or of its transpose in reverse mode, with the same solver, so the solver must suit both.

#type 'x t

The type for solvers of systems on values of type 'x.

#val dense : 'x t

dense materialises a as an n × n matrix, n the size of the vector, with one application to each vector of the standard basis, and solves it with Nx.solve: n + 1 applications with the check, and 2n³/3 flops. Its bound is c n ε (‖A‖ ‖u‖ + ‖r‖) + c n η (1 + ‖A‖): the backward error of an LU factorisation with partial pivoting, and the rounding of a solution of subnormal magnitude, whose components are multiples of η, with c = 16, ε the dtype's machine epsilon, η its least positive subnormal, ‖A‖ the Frobenius norm of the matrix and the vectors' norms Euclidean. It suits systems of up to some thousands of elements.

#val banded : width:int -> 'x t

banded ~width is for an a whose matrix has entries only within width of its diagonal, such as a finite difference in one dimension. It probes a with 2 width + 1 applications, each to the sum of the basis vectors of one residue modulo 2 width + 1, whose columns share no row, and factors the band by LU with partial pivoting, whose upper band grows to 2 width: 2 width + 2 applications with the check, and O(n width²) flops in loops of n steps. Its bound is dense's, with ‖A‖ the band's Frobenius norm. A band too narrow for a mixes its columns, and the solution misses the check.

Raises Invalid_argument if width < 0.

#val cg : rel:float -> budget:int -> precondition:('x -> 'x) -> 'x t

cg ~rel ~budget ~precondition is conjugate gradients from u = 0, preconditioned by precondition, for a symmetric positive-definite a such as a Hessian at a minimum or a damped Gauss–Newton matrix: an application of a and of precondition per iteration, and O(n) other work on three vectors. precondition approximates the inverse of a and is symmetric positive-definite itself; Fun.id is none. It stops when the residual it carries meets ‖a u − r‖ ≤ rel ‖r‖ + n η, the vectors' norms Euclidean and η the dtype's least positive subnormal, which is also its bound, so r = 0 gives u = 0 with the check alone. A solution of subnormal magnitude may miss it by its rounding, of order ‖a‖ n η, and end Stalled. budget iterations end the lane Budget_spent, and a direction p with pᵀ a p ≤ 0, which no positive-definite a has, ends it Stalled.

Raises Invalid_argument if rel is not in (0, 1) or if budget < 1.

val gmres : 
  restart:int ->
  rel:float ->
  budget:int ->
  precondition:('x -> 'x) ->
  'x t

gmres ~restart ~rel ~budget ~precondition is restarted GMRES from u = 0, preconditioned on the right by precondition, for any non-singular a. A cycle takes restart steps of Arnoldi's process, each an application of a and of precondition, orthogonalised twice by Gram–Schmidt, then the u of least residual in their span; it keeps restart + 1 vectors. Between cycles it applies a to u for the residual, and stops when ‖a u − r‖ ≤ rel ‖r‖ + n η, the vectors' norms Euclidean and η the dtype's least positive subnormal, which is also its bound (a solution of subnormal magnitude may miss it by its rounding, as cg's may), so r = 0 gives u = 0 with no application. After budget / restart cycles, the lane ends Budget_spent. precondition approximates the inverse of a and of its transpose, as a derivative solves both; Fun.id is none.

Raises Invalid_argument if restart < 1, if budget < restart or if rel is not in (0, 1).

val solve : 
  ('x, 'd) Nx.Ptree.t ->
  'x t ->
  ('x -> 'x) ->
  'x ->
  ('x, 'd) Solution.t

solve x s a r is the u with a u = r, found by s.

Raises Invalid_argument if the float tensors of r differ in dtype, or if a returns a value of another structure, dtype or shape than its argument.