GFT.jl — inverse generalized Fisher transformation

Official reference implementation for "Fast inversion of the generalized Fisher transformation of correlation matrices" (Archakov and Hansen; arXiv:2609.19028). A standard Julia package: standard libraries only (LinearAlgebra); Julia ≥ 1.6. The paper's numbers are from the release tagged v1.2.0.

using Pkg; Pkg.add(url = "https://github.com/reinhardhansen/GFT")
using GFT

Or run any of the scripts below directly from this folder; they activate the package environment themselves (julia -t 1 runtests.jl etc.). Pkg.test("GFT") runs the full test suite.

An R port with the same API (base R only, no dependencies) is on CRAN as GFT (install.packages("GFT")); its source lives in r/ (remotes::install_github("reinhardhansen/GFT", subdir = "r")). R package 1.2.0 ports this release function for function, including the variants, the quadrature preconditioner, inv_gft_path and gft_predict, inv_gft_anderson, and inv_gft_lbfgs.

• src/GFT.jl — the module: gft (forward map), inv_gft (GFT-FP+N, recommended: Newton steps for the log-diagonal residual, exact normalization along the vector of ones, quadratic forcing; the options residual = :gradient, forcing = :sqrt, normalize = false, phase = true, adaptive = false give the first-draft version, and exact_hess = true solves the Newton system with the explicit Hessian, the paper's full-Newton column; preconditioner = :quadrature or :auto replaces the diagonal preconditioner of the conjugate-gradient solve by the certified quadrature model of the paper's Section S6, :auto applying the selection rule), inv_gft_path (sequential inversion with warm starts and the tangent predictor gft_predict), and the comparison solvers inv_gft_fp (Archakov–Hansen fixed point), inv_gft_broyden (Chen–Fei–Yu; globalized = true adds the same initial fixed-point phase as inv_gft), inv_gft_newton (Newton in the form of Chen–Fei–Yu plus an Armijo line search with safeguard = false; the default adds the terminal safeguards of inv_gft), inv_gft_anderson (Anderson acceleration of the fixed point; guarded = true adds the quantified-decrease safeguard of the paper's Section 4, which makes it globally convergent) and inv_gft_lbfgs (limited-memory BFGS on the objective). All return an InvResult with fields x, C, iters, eighs, hvs, err, converged, hist.
• runtests.jl — test suite, including golden-value tests generated by an independent NumPy implementation (not distributed here), so a passing run is a cross-language verification.
• bench.jl — quick interactive benchmark (20 draws).
• overnight.jl — the full logged protocol behind the paper (1000 draws per random design, 500 at n=800, 1000 rate points); writes CSVs to results/.
• recheck_n800.jl — re-runs the 500 n=800 draws from the same RNG stream.
• final_checks.jl — same-input tolerance, same-initialization and safeguarded-Newton comparisons (supplement).
• check_serialized.jl — cross-language agreement on the serialized inputs in serialized/ (z vectors plus NumPy counts).
• figures.ipynb — Jupyter notebook (Julia kernel, needs Plots.jl) that builds the paper figure from results/, generates the LaTeX rows of Table 1, and prints every number quoted in the text.

using GFT
z = gft(C)          # C a nonsingular correlation matrix
r = inv_gft(z)      # r.C recovers C; r.converged, r.eighs, r.hvs

julia -t 1 runtests.jl        # verify first
julia -t 1 overnight.jl       # full logged protocol behind Table 1
julia -t 1 final_checks.jl    # supplement comparisons
julia -t 1 bench.jl extreme   # quick single blocks: toeplitz|random|extreme|big|warm

Notes for the paper's numbers: run single-threaded (-t 1; the script also sets BLAS.set_num_threads(1)), record hardware and Julia version, and ignore the first timing of any function (JIT compilation; the scripts warm up automatically). Random designs use MersenneTwister(18900217); draws are distributionally identical to, but not bit-identical with, the Python suite, so medians on random designs may differ slightly.

• v1.2.0 (September 2026): the revised GFT-FP+N of the JCGS submission. inv_gft solves H delta = -D ell (Newton for the log-diagonal residual) with forcing min(1/2, ||ell||), normalizes every evaluated point by the exact minimization of f along the vector of ones, drops the initial fixed-point phase, starts each line search from min(1, 2 t_prev), and tests descent before the untested-step rule (which now also requires ||ell||_inf < 1e-3). The first-draft algorithm remains available through keyword options and is the variant A of the ablation in ablation.jl (Section S5 of the paper). New: inv_gft_anderson(guarded = true), gft_predict, inv_gft_path, ablation.jl; overnight.jl adds a warm-start row with the predictor and final_checks.jl a guarded-Anderson column. Also new: the quadrature preconditioner (preconditioner option of inv_gft, with quadrature_precond, logphi, choose_order) and precond.jl, which produces the paper's Tables S5 and S6 (results/precond.csv, precond_rule.csv, precond_warm.csv). serialized/py_counts.csv holds the NumPy counts of the revised algorithm. Median eigendecompositions fall from 17/23 to 10/14 on the z-designs and by two or three on the structured designs, with no failures on any design including n = 800.
• v1.1.0 (September 18, 2026): inv_gft gains exact_hess, inv_gft_newton gains safeguard (default true; false is the published comparator), new comparators inv_gft_anderson and inv_gft_lbfgs; overnight.jl runs the full-Newton column (newtonsafe), final_checks.jl records wall times and runs the Anderson/L-BFGS comparison, recheck_n800.jl takes the tolerance as an argument. results/ holds the logged rerun of the complete protocol on an Apple M2 Max under Julia 1.13.0 (results/log.txt, results/machine.txt, results/*_run.log); every GFT-FP, Broyden and GFT-FP+N count of the v1.0.1 run is reproduced exactly or within one eigendecomposition, the counts of the published Newton variant moved at the rounding floor as documented in the paper, and one n = 800 draw fails at tolerance 1e-13 because the stopping criterion cannot be evaluated below 2e-13 on it (all 500 converge at 1e-12, results/n800_tol1.0e-12.csv).
• v1.0.1 (July 27, 2026): the full-Newton comparator's default iteration cap is corrected from an unintended 100 to the 500 stated in the paper's protocol, and results/ holds the logged rerun under the corrected cap. The rerun reproduced every iteration count and failure total of the v1.0.0 run exactly; only wall times differ, within run-to-run noise.
• v1.0.0 (July 19, 2026): initial release. Its logged results were produced with the full-Newton comparator capped at 100 iterations; all other methods used the documented caps.

Known sensitivity: the iteration count of Newton with Armijo backtracking (inv_gft_newton with safeguard = false) on near-singular deterministic designs (e.g. Toeplitz rho=0.99) is unstable to roundoff-level perturbations — accept/reject decisions of its Armijo line search sit near boundaries and flip across BLAS builds or algebraically equivalent reformulations (observed counts 7, 10, 23 and 34 for the same problem). The counts of GFT-FP, Broyden, GFT-FP+N and the safeguarded solvers are stable. This fragility is intrinsic to Newton with an objective-based line search at tight tolerance without the untested-step rule, and is discussed in the paper.