EvoCG: Self-Evolving Precision Control for Faster Conjugate Gradient

arXiv:2611.04469v1 [cs.DC] 3 Nov 2026

EvoCG: Self-Evolving Precision Control for Faster Conjugate Gradient Jiaqi Zhai∗†‡

Haozhi Han∗†

Kaiyi Huang‡

jqzhai@hust.edu.cn Huazhong University of Science and Technology Wuhan, China

haozhi.han@stu.pku.edu.cn School of Computer Science, Peking University Beijing, China

huangky@hust.edu.cn Huazhong University of Science and Technology Wuhan, China

Wenju Zhao‡

Yifeng Chen

Yunquan Zhang

wjzh@hust.edu.cn Huazhong University of Science and Technology Wuhan, China

cyf@pku.edu.cn School of Computer Science, Peking University Beijing, China

zyq@ict.ac.cn Institute of Computing Technology, Chinese Academy of Sciences Beijing, China

Ting Cao

Yunxin Liu

Hai Jin‡

tingcao@mail.tsinghua.edu.cn Institute for AI Industry Research, Tsinghua University Beijing, China

liuyunxin@air.tsinghua.edu.cn Institute for AI Industry Research, Tsinghua University Beijing, China

hjin@hust.edu.cn Huazhong University of Science and Technology Wuhan, China

Xuanhua Shi✉‡

Kun Li✉

xhshi@hust.edu.cn Huazhong University of Science and Technology Wuhan, China

likungw@gmail.com Institute for AI Industry Research, Tsinghua University Beijing, China

Abstract Modern GPUs offer far higher throughput in low precision than in FP64, but exploiting this in conjugate gradient (CG) solvers is difficult because tolerance to reduced precision evolves with the Krylov state. We present EvoCG, a mixed-precision CG solver that treats tile-wise precision as an online control action rather than a fixed setting: a reinforcement-learning policy in the iteration loop selects per-tile precision from the evolving solver state, made profitable by co-design with a Tile-CSC unified mixed-precision SpMV kernel and a low-overhead in-loop runtime. Across 47 SPD matrices from eight scientific domains, EvoCG is the fastest solver overall, with geometric-mean speedups of 2.52×, 3.06×, 2.48×, 1.93×, and 6.99× over cuSPARSE, PETSc, Ginkgo, Mille-FP64, and MilleMixed. On the same substrate, it improves time-to-solution over EvoCG-FP64 by 1.28× at only 11.4% more iterations and 1.09% overhead, with FP64-recomputed final true-residual norms nearly identical to those of FP64, and extends to Jacobi-preconditioned ∗ Both authors contributed equally to this research. † Work done during an internship at Institute for AI Industry Research(AIR), Tsinghua

University. ‡ National Engineering Research Center for Big Data Technology and System; Services Computing Technology and System Lab; Cluster and Grid Computing Lab; School of Computer Science and Technology, Huazhong University of Science and Technology, Wuhan 430074, China. PPoPP’27, Salt Lake City, Utah 2026. ACM ISBN 978-x-xxxx-xxxx-x/YYYY/MM https://doi.org/10.1145/nnnnnnn.nnnnnnn

CG without retraining, yielding a further 1.15× geometric-mean speedup over EvoCG-FP64.

Keywords Conjugate Gradient, Mixed-Precision Computing, Reinforcement Learning, GPU Acceleration

1

Introduction

The rise of AI has reshaped the precision landscape of modern processors [23, 29]: contemporary GPUs deliver only tens of TFLOPS in FP64, whereas lower-precision formats such as FP32, TF32, FP16, BF16, and INT8 offer peak throughputs one to two orders of magnitude higher [15, 16]. Scientific computing, in contrast, still relies predominantly on uniform double precision for robustness and convergence guarantees [19], creating a fundamental efficiency gap: massive low-precision capacity that traditional HPC workloads cannot exploit safely and systematically [42]. This mismatch is pronounced in Krylov subspace solvers [20, 24, 33]. Among them, the Conjugate Gradient (CG) method is a cornerstone HPC sparse solver underpinning finite element analysis [11], computational fluid dynamics [13], and computational electromagnetics [41]. Yet production CG implementations—including cuSPARSE-based solvers [31], PETSc [5], and Ginkgo [3]—still execute the entire iteration uniformly in FP64, leaving the abundant low-precision throughput of modern GPUs untapped.

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

Within the CG iteration, Sparse Matrix–Vector Multiplication (SpMV) dominates runtime, often 60%–95% of time-to-solution for large sparse systems [6, 27, 43]. Lowering SpMV precision is thus a major lever for accelerating CG. Our experiments show that, on average, 56.9% of SpMV computation across the 47 matrices is executed in FP32 or BF16. The challenge is that reduced-precision SpMV perturbs the Krylov recurrence [7, 35], and CG’s tolerance to such perturbations evolves with the solver state, potentially slowing or even destabilizing convergence [8, 12, 25, 46]. Mixed-precision CG has been explored in several forms, yet precision remains largely an open-loop decision: early approaches confine low precision to preconditioners or inner solvers, keeping SpMV in high precision [26, 38]. Recent GPU-oriented methods lower SpMV precision but follow static or rule-based schedules agnostic to the Krylov state [46]—schedules that, in our evaluation, frequently fail to converge on sensitive problems (Sec. 4.7). RL-based approaches select precision as coarse-grained offline autotuning, with no execution substrate to turn low precision into actual speedup [12]. Existing methods thus expose precision heterogeneity without closing the feedback loop or delivering its performance benefit. These limitations point to a deeper question: Can Krylov solvers become self-evolving systems in which precision is adaptively controlled throughout the iteration process? We answer this question with EvoCG, built on a key insight: precision is control, not configuration—not fixed offline as a static setting, but regulated online as a solver action driven by the evolving Krylov state. EvoCG realizes this as a closed loop: a reinforcement-learning policy embedded in the CG iteration reads the Krylov state and selects the tile-wise precision, whose induced inexactness in turn shapes the next state. Precision thus becomes a controllable input, and inexactness a managed degree of freedom that balances convergence and hardware efficiency. Turning low precision into end-to-end solver speedup requires meeting three requirements simultaneously: the solver must recognize when reduced precision is tolerable for the current Krylov state, convert that tolerance into real speedup despite tile-level precision heterogeneity, and integrate this control at overhead low enough for an iterative solver—missing any one forfeits the benefit. EvoCG addresses these with three co-designed modules: EvoCtrl maps the current search direction 𝑝𝑘 to a tile-level numerical state and selects a tile-wise precision assignment over {FP64, FP32, BF16} with a learned policy; EvoSpMV executes the resulting heterogeneous SpMV in a single warp-centric kernel built on Tile-CSC, a storage format co-designed to absorb tile-level precision heterogeneity without control-flow divergence; and EvoRun amortizes the control overhead on the critical path through device-side inference, interval-based reuse, and CUDA-Graph replay. On 47 SPD SuiteSparse matrices spanning eight scientific domains, EvoCG is the fastest solver overall, with geometric-mean speedups of 2.52×, 3.06×, 2.48×, 1.93×, and 6.99× over cuSPARSE, PETSc, Ginkgo, Mille-FP64, and Mille-Mixed. Against EvoCG-FP64, our FP64 CG solver on the identical substrate, dynamic precision control improves end-to-end performance by 1.28× on average at only 11.4% more iterations, and the average policy overhead is only 1.09%. The same policy further extends to Jacobi-preconditioned

Zhai et al.

CG without retraining, delivering a 1.15× geometric-mean speedup over EvoCG-FP64. This paper makes the following contributions: • We exploit CG’s state-dependent precision tolerance to formulate mixed-precision CG as a reinforcement-learning control problem, turning tile-wise precision into a learned solver action that yields end-to-end speedup. • We turn tile-wise low-precision decisions into solver-level speedup with Tile-CSC and a unified warp-centric SpMV kernel that runs mixed-precision tiles without control-flow divergence. • We keep precision control on the CG critical path at low overhead, via device-side inference, interval-based reuse, and CUDA-Graph-replayed fused execution.

2 Background and Challenges 2.1 Conjugate Gradient The Conjugate Gradient (CG) method solves symmetric positive definite linear systems 𝐴𝑥 = 𝑏 by iteratively constructing Krylov subspaces. At iteration 𝑘, CG maintains a solver state 𝑧𝑘 = {𝑥𝑘 , 𝑟𝑘 , 𝑝𝑘 },

(1)

where 𝑥𝑘 is the approximate solution, 𝑟𝑘 the residual, and 𝑝𝑘 the search direction. CG evolves this state through the short recurrence shown in Fig. 1, consisting of one SpMV, two dot products, and three AXPY operations per iteration. Among these operations, SpMV dominates execution time: Fig. 2 shows that across 47 sparse systems from eight domains, it accounts for 67.1% of runtime on average and exceeds 80% for many large problems. Reducing SpMV precision improves arithmetic efficiency but introduces numerical perturbations into the CG recurrence: 𝑦˜𝑘 = (𝐴 + 𝐸𝑘 )𝑝𝑘 ,

(2)

where 𝐸𝑘 represents the error introduced by reduced-precision arithmetic. Inexact Krylov theory interprets this as applying CG to a sequence of perturbed operators and shows that such perturbations affect both attainable accuracy and convergence behavior [7, 35, 40]. The impact on attainable accuracy is characterized by the residual gap between the true residual 𝑟𝑚 = 𝑏 − 𝐴𝑥˜𝑚 and the recursively updated residual 𝑟˜𝑚 : ∥𝑟𝑚 − 𝑟˜𝑚 ∥ ≤

𝑚 ∑︁

∥𝐸𝑘 ∥|𝜂𝑘(𝑚) |,

(3)

𝑘=1

SpMV

Search Direction

High Computation Cost

Solution

Residual

Figure 1: One iteration of the CG recurrence. SpMV is the dominant cost.

EvoCG : Self-Evolving Precision Control for Faster Conjugate Gradient

Figure 2: Dominance of SpMV in CG execution across 47 sparse linear systems from eight scientific domains. where 𝜂𝑘(𝑚) represents the contribution of the 𝑘-th perturbation to the residual gap [35]. Therefore, controlled perturbations can preserve the attainable final accuracy. For short-recurrence methods such as CG, perturbations may further alter the convergence trajectory through residual drift and loss of orthogonality [21, 25, 40], introducing a fundamental tradeoff: lower precision reduces the cost of each SpMV, but the added perturbations may increase the total iteration count. Inexact Krylov theory characterizes the relationship between perturbations and convergence, but it does not determine which precision strategy minimizes end-to-end solution time. This motivates adaptive precision control, where each precision decision must balance immediate execution cost and future convergence behavior.

2.2

Reinforcement Learning

The above precision-control problem naturally forms a sequential decision process. At each CG iteration, a controller observes the evolving solver state, selects precision actions, and receives feedback from both execution cost and convergence behavior. Unlike static precision assignment or residual-only heuristics, the optimal decision depends on long-term effects across iterations. Reinforcement learning provides a natural framework for this setting by optimizing policies through long-term rewards. Proximal Policy Optimization (PPO) [34] is adopted because its clipped policy updates improve optimization stability for long-horizon control problems [44].

2.3

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

C1: Capturing long-term numerical effects of precision decisions. Precision choices affect not only current SpMV cost but also future convergence through the Krylov recurrence, making it difficult to design states and rewards that capture the accuracyperformance tradeoff. C2: Performance-effective execution under sparsity and heterogeneous precision. Fine-grained precision control introduces irregular computation, whereas GPUs favor regular execution, requiring cross-layer co-design of sparse format, kernel mapping, and scheduling. C3: Low-overhead integration of adaptive control. A learned precision policy must be integrated into the CG execution loop with minimal runtime overhead, requiring efficient inference and solveraware execution mechanisms.

3 EvoCG 3.1 Overview As shown in Fig. 3, EvoCG comprises three tightly coupled modules: EvoCtrl (control) learns an RL policy that selects tile-wise precision from the evolving Krylov state. EvoSpMV (execution) turns these decisions into an efficient sparse substrate via the co-design of Tile-CSC and a unified warp-centric mixed-precision kernel. EvoRun (runtime) embeds policy execution into the CG loop via low-overhead inference, interval reuse, and fused execution. Together, they determine what precision to use, how to execute it efficiently, and when this control enters the solver.

3.2

EvoCtrl: Closed-Loop Precision Control

3.2.1 Control Objective Formulation. EvoCtrl minimizes CG timeto-solution while preserving convergence robustness. Let 𝑇 ({𝑎𝑘 }) be the total execution time to reach a prescribed residual tolerance under precision decisions {𝑎𝑘 }: 𝑇 ({𝑎𝑘 }) =

𝐾 ( {𝑎 𝑘 } ) −1 ∑︁

𝐶 (𝑎𝑘 ),

where 𝑎𝑘 is the precision configuration at iteration 𝑘, 𝐶 (𝑎𝑘 ) is its per-iteration cost, and 𝐾 ({𝑎𝑘 }) is the iteration count. EvoCtrl

Offline training

Online Interval Update

EvoSpMV

Precision Control Formulation PPO Quantization BF16

Solver State

column-block

CG maps naturally onto an MDP. The evolving solver state {𝑥𝑘 , 𝑟𝑘 , 𝑝𝑘 } supplies a Markov state, tile-wise precision forms a discrete action space, and the CG recurrence itself is the transition dynamics, requiring no learned environment model. The one non-obvious element is the reward, which must encode a delayed objective: a cheap low-precision SpMV now may cost extra iterations later, so end-to-end time-to-solution, not per-iteration cost, is what the policy must optimize. Sec. 3.2 instantiates each component.

EvoRun Low-overhead inference

Tile-CSC

Interval Reuse Scheduling FP64 Fused CG update Compute Update Compute Update

2.4

Challenges

Unlocking low precision in CG while preserving convergence robustness poses three challenges.

(4)

𝑘=0

Figure 3: EvoCG overview.

Mix-Precision SpMV Kernel

Speedup Error Injection

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

Zhai et al.

Making this dependence explicit, mixed-precision SpMV is modeled as: 𝑦˜𝑘 = (𝐴 + 𝐸𝑘 (𝑎𝑘 ))𝑝𝑘 , (5) where 𝐸𝑘 (𝑎𝑘 ) denotes the perturbation induced by action 𝑎𝑘 . This yields the controlled solver dynamics: 𝑧𝑘+1 = 𝐹 (𝑧𝑘 , 𝑎𝑘 ),

(6)

where 𝑧𝑘 denotes the solver state. Precision decisions therefore affect both the immediate execution cost and the future convergence trajectory, making EvoCtrl recurrence-driven control rather than a static configuration rule. 3.2.2 Precision Control MDP. EvoCG models the controlled CG process as an MDP (S, A, P, R, 𝛾): each CG iteration is a decision step and a complete solve an episode. State Space S. The MDP state is constructed from the CG search direction 𝑝𝑘 , motivated by two properties. First, through the CG recurrence, 𝑝𝑘+1 = 𝑟𝑘+1 + 𝛽𝑘 𝑝𝑘 , (7) 𝑝𝑘 carries accumulated information about the convergence regime. Second, 𝑝𝑘 is the direct input to SpMV, so its tile-level distribution determines the magnitude of partial sums and their sensitivity to reduced precision. These two properties make 𝑝𝑘 a well-suited state for convergence-aware control. To expose both global solver dynamics and localized computation, 𝑝𝑘 ∈ R𝑛 is partitioned into 𝑛 tiles segments of size 𝐶, each normalized to form the policy input:

INIT: policy setup and FP64 baseline estimation. ENV: RL controller operations. SOVR: CG recurrence updates. SpMV: tile-wise mixed-precision SpMV execution.

ENV Input: Matrix 𝐴, right-hand side 𝑏, tile size 𝐶 ENV Output: Trained policy network 𝜋𝜃 3 INIT Initialize policy 𝜋𝜃 , value network 𝑉𝜃 ; 4 INIT Run full FP64 CG to obtain 𝑏𝑎𝑠𝑒_𝑖𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛; 5 ENV 𝑚𝑎𝑥_𝑖𝑡𝑒𝑟 ← 1.5 × 𝑏𝑎𝑠𝑒_𝑖𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛; 6 ENV 𝑛𝑡𝑖𝑙𝑒𝑠 ← ( ⌈𝐴.shape[1]/𝐶 ⌉ ); 7 for 𝑒𝑝𝑖𝑠𝑜𝑑𝑒 = 1 to 1000 do 8 SOVR Initialize 𝑥 ← 0, 𝑟 ← 𝑏 − 𝐴𝑥, 𝑝 ← 𝑟 ; 9 for 𝑘 = 0 until ∥𝑟 ∥/∥𝑟 0 ∥ ≤ 𝜏 or 𝑘 = max_iter do 10 for 𝑖 ← 0 to 𝑛 tiles − 1 do 11 ENV 𝑝 (𝑖 ) ← 𝑝 [𝑖𝐶 : (𝑖 + 1)𝐶 − 1]; 12 ENV 𝑠 (𝑖 ) ← 𝜙 (𝑝 (𝑖 ) ); 13 ENV 𝑎 (𝑖 ) ∼ 𝜋𝜃 (𝑠 (𝑖 ) ); 1 2

15

for each tile 𝑖 do SpMV 𝑦 ← 𝑦 + 𝐴 (:,𝑖 ) 𝑝 (𝑖 )

16

SOVR 𝛼 ← (𝑟 𝑇 𝑟 )/(𝑝𝑇 𝑦);

17

20

SOVR 𝑥 ← 𝑥 + 𝛼𝑝; SOVR 𝑟 new ← 𝑟 − 𝛼 𝑦; 𝑇 𝑟 𝑇 SOVR 𝛽 ← (𝑟 new new )/(𝑟 𝑟 ); SOVR 𝑝 ← 𝑟 new + 𝛽𝑝;

21

ENV

22

SOVR 𝑟 ← 𝑟 new ;

14

18 19

23

𝑝𝑘(𝑖 ) = 𝑝𝑘 [𝑖𝐶 : (𝑖 + 1)𝐶 − 1], 𝜙 (𝑝𝑘(𝑖 ) ) = 𝑝𝑘(𝑖 ) /max(∥𝑝𝑘(𝑖 ) ∥ 2, 𝛿),

Algorithm 1: EvoCG Training Phase.

24

ENV ENV

(by 𝑎 (𝑖 ) );

  (𝑎) ∥ 𝑅 ← 𝑤1 ln ∥𝑟∥𝑟 − 𝑤2 𝐶 + 𝑤3 ⊮ ∥𝑟∥𝑟new∥ ∥ ≤ 𝜏 ; 𝑛 new ∥ tiles

0

Update 𝜋𝜃 and 𝑉𝜃 using PPO objective;

return 𝜋𝜃

(8) from stability-oriented defaults and reused unchanged across all 47 matrices, with the convergence terms weighted higher because early low-precision exploration can destabilize CG.

(𝑛 −1) 𝑠𝑘 = (𝜙 (𝑝𝑘(0) ), . . . , 𝜙 (𝑝𝑘 tiles )),

where 𝛿 prevents division by zero. Action Space A. At iteration 𝑘, EvoCtrl outputs a tile-wise precision assignment vector: (𝑛

𝑎𝑘 = (𝑎𝑘(0) , . . . , 𝑎𝑘 tiles

−1)

),

(9)

where 𝑎𝑘(𝑖 ) ∈ {0, 1, 2} encodes the precision {FP64, FP32, BF16} for the vector tile 𝑝𝑘(𝑖 ) and its associated column block, as defined in Eq. (10): 𝐴 (:,𝑖 ) = 𝐴[:, 𝑖𝐶 : (𝑖 + 1)𝐶 − 1],

(10)

Reward Function R. The reward is designed as an objective for minimizing end-to-end time-to-solution. It jointly captures convergence progress, execution cost, and successful termination:   ∥𝑟𝑘 −1 ∥ 𝐶 (𝑎𝑘 ) ∥𝑟𝑘 ∥ 𝑅𝑘 = 𝑤 1 ln − 𝑤2 + 𝑤3⊮ ≤𝜏 , (11) ∥𝑟𝑘 ∥ 𝑛 tiles ∥𝑟 0 ∥ where 𝜏 is the convergence tolerance. The first term rewards residual reduction, its logarithmic form normalizing improvement across iterations. The second term penalizes per-iteration cost 𝐶 (𝑎𝑘 ) through a scale-independent per-tile signal, letting the same reward generalize across matrix scales. The terminal bonus prevents aggressive low-precision choices that cut immediate cost but miss the target accuracy. We set 𝑤 1 = 10, 𝑤 2 = 0.2, and 𝑤 3 = 30, chosen

Transition Dynamics P and Discount Factor 𝛾. The transition dynamics are the CG recurrence itself, so PPO interacts directly with the mixed-precision solver as the environment without an explicit model of P. We set 𝛾 = 0.99.

3.2.3 Policy Training. EvoCtrl is instantiated as a lightweight MLP parameterizing both policy and value: each normalized segment 𝜙 (𝑝𝑘(𝑖 ) ) ∈ R𝐶 passes through a two-layer ReLU backbone. The policy head outputs a three-class distribution over {FP64, FP32, BF16} yielding 𝑎𝑘(𝑖 ) , and the value head predicts the expected return. The network is trained with PPO in a solver-in-the-loop manner, summarized in Algorithm 1: lines 10–13 construct tile states and sample tile-wise actions, line 15 executes the mixed-precision SpMV inside the CG recurrence, the SOVR lines perform the standard CG updates, line 21 evaluates the reward of Eq. (11), and line 23 stores solver transitions and apply the clipped PPO update whenever 1024 transitions have been collected. The tile-wise precision actions learned by EvoCtrl are later consumed by EvoSpMV, which serves as the sparse execution substrate for efficient GPU realization.

EvoCG : Self-Evolving Precision Control for Faster Conjugate Gradient

3.3

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

EvoSpMV: Mixed-Precision Sparse Execution

EvoSpMV is the execution module of EvoCG. It executes EvoCtrl’s tile-wise precision decisions under irregular sparsity without fragmenting into precision-specific kernel paths. The challenge is that mixed precision makes SpMV heterogeneous in both sparsity and arithmetic cost. EvoSpMV addresses this by co-designing the Tile-CSC layout with a unified warp-centric mixed-precision kernel, confining precision heterogeneity to tiles while keeping execution regular at the warp-group and lane levels.

BF16

FP64

FP32

warps group2 coarse-grained balance

BF16

BF16

FP32

FP64

3 N/R

7 ... N/C

BF16

2

FP32 FP32

FP32

1

BF16

0

warps group1 warps group4 3 4 5 6

BF16

1 2 FP64 FP64

0

CSR inter-tile tileRowPtr: [0, 3, 6, 9, 12, 15] tileColIdx: [1, 7, N/C, 1, 3, 5, 2, 4, 8, 4, 6, N/C, 0, 3, 7] tileInfo: [CSC intra-tile meta data] scale: [ ]

(a) Inter-Tile Row-Block Scheduling

lane0

lane1

lane31 Memory Coalescing

0 1 2 3 4 5 ...

...

...

R-1

R

padding 0 1 ... 31 32 33 ... 63 colPtr: rowIdx: FP64-val: FP32-val: BF16-val:

0 1 ... 31 32 33 ... 63

lane-level balance lane 0

lanes

CSC intra-tile

warps group0 warps group3

lane 1

...

... lane 31

[0, 4, 6, ... , x+2] sentinel row index [[0, 3, R-1, R], [2, 5], ... , [2, 4]]

[[ , , , ], [ , ,], ... , [ , ]] [[ , , , ], [ , ,], ... , [ , ]] [[ , , , ], [ , ,], ... , [ , ]]

(b) Intra-Tile CSC Execution

Figure 4: EvoSpMV: co-design of Tile-CSC and a warp-centric mixed-precision SpMV kernel. As shown in Fig. 4(a), the sparse matrix is partitioned into 𝑅 × 𝐶 tiles, where 𝐶 is fixed to 64 and 𝑅 is determined during preprocessing. Each tile corresponds to the intersection between an 𝑅-row block and a consecutive 64-column block, serving as the unified granularity for Tile-CSC storage, EvoSpMV scheduling, and runtime precision control. 3.3.1 Inter-Tile Row-Block Scheduling. The inter-tile challenge is balancing warp-group workload under heterogeneous tile precisions. EvoCG organizes non-empty tiles row-block-wise in CSR style, where tileRowPtr delimits each row-block’s tile interval and tileColIdx records column-block indices (Fig. 4(a)). Each rowblock is assigned a fixed group of 𝑊 warps traversing its non-empty tiles in a strided manner (warp 𝑤 processes tiles 𝑡 0 + 𝑤, 𝑡 0 + 𝑤 + 𝑊 , . . .). This CSR-style organization disperses tiles of different precision across warp groups: each group spans a mixture of FP64, FP32, and BF16 tiles rather than a precision-homogeneous stream, averaging the precision-dependent cost variation into a uniform workload. A CSC-style layout would instead cluster same-precision tiles within a warp group, causing significant inter-warp imbalance. 3.3.2 Intra-Tile CSC Execution. The intra-tile challenge is dynamic precision selection without sacrificing warp-level regularity. As Fig. 4(b) shows, EvoSpMV encodes each tile in CSC format and assigns one warp per tile. The tile descriptor tileInfo records the global-memory offsets of colPtr, rowIdx, and the value streams: within a tile, colPtr (length 𝐶 + 1) defines each column’s nnz range, and rowIdx stores row offsets relative to the row-block base. The warp traverses the 64 columns cooperatively, lane 𝑖 processing columns 𝑖 and 𝑖 + 32. Accesses to x, scale, and colPtr coalesce

across lanes, while rowIdx and value accesses stay column-local, confining sparse irregularity to column-local traversal rather than warp-wide divergence. This column-centric path aligns naturally with the mixed-precision layout. Quantization is applied per column, so all nonzeros in a column share one scaling factor, and column traversal amortizes the loads of 𝑥 [ 𝑗] and 𝑠 [ 𝑗]. Since FP32 and BF16 share dynamic range, one column-wise factor serves both. To switch precision without altering traversal, each nonzero carries three value streams (FP64-val, FP32-val, BF16-val) over the same structural indices, so precision selection changes only which stream is read. Columns are padded to a multiple of two with sentinel-marked zeros, preserving vectorized access and structural alignment without affecting logical sparsity. EvoSpMV thus runs heterogeneous tile precisions in one indexing framework with no kernel-switch overhead, providing the stable backbone EvoRun embeds into the CG critical path. 3.3.3 Cross-Level Balance Tuning. The row-block height 𝑅 is EvoSpMV’s main preprocessing knob for cross-level balance. EvoCG tunes 𝑅 over {32, 64, . . . , 16384} to even non-empty tile counts across rowblocks for inter-tile warp-group balance, which in turn reshapes the nnz distribution within 64-column tiles for lane-level balance.

3.4

EvoRun: Low-Overhead In-Loop Runtime

As illustrated in Fig. 5, EvoRun is the runtime module of EvoCG. It integrates EvoCtrl and EvoSpMV directly into the CG iteration loop. Policy inference remains part of the critical path, but its cost is reduced and amortized through a fused device-side inference kernel and interval-based action reuse. The remaining CG operators are fused and replayed through a CUDA Graph, reducing launch and intermediate-memory overhead around the controlled SpMV. Policy Inference. Since policy execution runs inside CG iteration, its cost must stay minimal. We fuse the whole policy inference into one kernel, avoiding intermediate write-backs and extra launches. Policy weights are quantized offline to FP16 and reordered into a WMMA-friendly layout, enabling FP16 Tensor Core inference with FP32 accumulation. Interval Reuse Scheduling. Two observations motivate reusing precision decisions across iterations. First, column-block precision patterns exhibit strong temporal persistence between adjacent iterations, so per-iteration policy invocation is largely redundant. Second, precision transitions are asymmetric: delaying an upward transition may leave a block insufficiently accurate and even break convergence, whereas delaying a downward one only costs efficiency. Together, these suggest that the policy can be invoked periodically, provided upward transitions are not unduly delayed. Guided by this, we set the reuse interval 𝐾 ★ offline. From a periteration run of EvoCG, we record each column block’s precision and, over candidate intervals 1, 2, 4, 8, . . . , measure the fraction of upward transitions that reuse would delay. 𝐾 ★ is the largest interval keeping this fraction below 0.1. Online, the policy runs every 𝐾 ★ iterations, reusing its decision within each window.

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

Zhai et al.

...

Offline

criterion for K*

FP32 BF16

Online Inference Precision Control applied for K* iterations

Interval Reuse Scheduling i=0

Optimal

i=K*

CG iteration

Profile-guided Asymmetric Transition Management (Preserve upward transitions)

convergence time

ReLU

EvoSpMV (Block

)

Fused CG Update Kernel Phase 0

Phase 1

Vectorized Data Load

Compute

Reduce

Phase 2

Phase 3

Update Write

FMA Reduce

Compute & Update Write Grid-wide Sync 2

...

CUDA Graph Boundary

Grid-wide Sync 1

FP64 ...

Shared Mem Global Mem

FP16 weights are executed on TensorCore

Register

Offline Quantization

EvoCtrl: Policy Network

column-block

Precision Control tilesize=64

FMA

i=2K* ...

K*

Figure 5: EvoRun: system-level realization of EvoCG with in-loop policy inference—GPU-resident policy execution, action reuse over 𝐾 ★ iterations, and CUDA-Graph replay of the fused solver backbone. Fused Execution and Graph Replay. To complete the runtime path, EvoRun compresses the remaining CG operators—dot-product reductions, residual update, and search-direction update—into a single cooperative fused-step kernel that completes global reduction, scalar broadcast, and vector updates entirely on device via grid-wide synchronization, removing repeated launches, redundant memory traffic, and intermediate write-backs. The SpMV + fused-step backbone is further captured in a CUDA Graph to eliminate multikernel launch overhead, while policy inference stays outside the graph and triggers only once every 𝐾 ★ iterations.

4 Experiments 4.1 Experimental Setup Hardware and software. All experiments run on an NVIDIA A100 (40 GB HBM) with an AMD EPYC 7K62 CPU (8 cores, 16 threads) and 93 GiB host memory, under Ubuntu 22.04, CUDA 12.2 (driver 535.216.01), GCC 11.4.0, CMake 3.22.1. Baselines. We compare EvoCG against five baseline configurations from four solver systems: PETSc v3.24.5 [30] , Ginkgo v1.11.0 [3], an FP64 CG solver based on cuSPARSE [31] and cuBLAS, and two configurations of Mille-feuille [46]. The cuSPARSE baseline uses CSR SpMV from cuSPARSE-cu12 v12.3.1 together with cuBLAS for the remaining vector operations. We denote the two Mille-feuille configurations as Mille-FP64 and Mille-Mixed. MilleFP64 executes the entire solver in FP64, whereas Mille-Mixed uses the released mixed-precision implementation. PETSc, Ginkgo, and the cuSPARSE baseline do not support mixed-precision CG in our experimental setting and therefore run entirely in FP64. Precision configuration. The tile-level action space includes FP64, FP32, and BF16. FP8 is excluded because A100 offers no practically beneficial FP8 path for our workflow, and FP8 actions destabilize RL training and prevent reliable policy convergence.

Unified solver configuration. All methods use the same configuration: initial guess 𝑥 0 = 0, right-hand side 𝑏 = 𝐴1 (exact solution 1), the recursively updated relative-residual stopping criterion ∥𝑟𝑘 ∥ 2 /∥𝑟 0 ∥ 2 ≤ 10−10 , and a maximum of 10,000 iterations. Datasets. We evaluate EvoCG on 47 sparse matrices from the SuiteSparse Matrix Collection [18], covering eight scientific domains. Our benchmark set includes all SuiteSparse matrices that meet three criteria: they are symmetric positive definite (SPD), have dimension greater than 10,000, and converge within 10,000 CG iterations. The dimension threshold ensures GPU utilization and SpMV dominance. The iteration threshold keeps solver-in-the-loop training feasible. The set was not filtered by conditioning and includes difficult instances such as nd6k and cbuckle (estimated condition numbers > 107 ). Training configuration. The precision-control policy is trained with PPO using 𝛾 = 0.99, learning rate 10−4 , 10 epochs per update, minibatch size 128, clip ratio 0.2, value-loss coefficient 0.25, entropy coefficient 0.02, update interval 1024, hidden sizes [64, 64], gradient clipping 0.5, advantage standardization, and GAE parameter 𝜆 = 0.95. These hyperparameters follow stability-oriented defaults [34] refined by a single small validation sweep, shared unchanged across all matrices and never tuned per matrix. Training protocol. EvoCG-Mixed uses a matrix-specific policy trained offline on the target matrix. This setting targets repeatedsolve deployments, where training is a one-time cost amortized over subsequent solves. Sec. 4.5 evaluates transferring one mesh2e1trained policy to unseen matrices. Reported times cover the online solve only, excluding offline policy training, whose cost is discussed in Sec. 4.8.

4.2

End-to-End Solver Performance

As shown in Fig. 6, under one common configuration and convergence criterion, EvoCG is consistently the fastest solver across the 47 SPD matrices, with geometric-mean speedups of 2.52×, 3.06×,

EvoCG : Self-Evolving Precision Control for Faster Conjugate Gradient

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

Figure 6: End-to-end performance vs. cuSPARSE, PETSc, Ginkgo, and two Mille-feuille variants across 47 SPD matrices: achieved GFLOPS (top) and solver-level speedup of EvoCG over each baseline on a log10 axis (bottom), against nonzeros (log10 ). 2.48×, 1.93×, and 6.99× over cuSPARSE, PETSc, Ginkgo, MilleFP64, and Mille-Mixed, respectively. The corresponding maximum speedups reach 4.47×, 4.89×, 4.44×, 4.87×, and 22.49×. For MilleFP64 and Mille-Mixed, we report results only on matrices where their released implementations complete successfully. Against cuSPARSE, PETSc, and Ginkgo, the gain comes from full-stack co-design: adaptive mixed precision cuts arithmetic cost whenever the Krylov state tolerates inexactness, Tile-CSC and the warp-centric kernel absorb the execution irregularity of sparsity plus precision heterogeneity, and fused execution removes redundant launches and intermediate memory traffic. Crucially, none of these gains is erased by convergence degradation. The Mille-FP64 comparison is particularly informative: MilleFP64 often outperforms cuSPARSE, PETSc, and Ginkgo thanks to its fused-kernel design, yet it remains a fixed-precision FP64 solver without EvoCG’s co-design around precision heterogeneity, so EvoCG achieves a clear 1.93× geometric-mean speedup. EvoCG achieves an even larger speedup over Mille-Mixed, and the reasons are instructive. The mega-kernel design of Mille-Mixed introduces device-side spin-wait synchronization, four precision copies, periteration casts, and d_vis-guided warp branches. EvoCG avoids these overheads with a staged pipeline, in which Tile-CSC mixedprecision SpMV is followed by a cooperative fused kernel synchronized grid-wide via cooperative_groups. This contrast illustrates a broader point: mixed precision brings not only lower arithmetic cost but also execution heterogeneity, since tiles of different precision incur different amounts of work. Mixed precision thus pays off only when the storage format, kernel mapping, and solver execution are co-designed to absorb this imbalance.

4.3

Ablation of Dynamic Precision Control

We next isolate dynamic precision control. The reference EvoCGFP64 is our CG solver on the identical Tile-CSC substrate, fused kernels, and CUDA-Graph runtime, with every tile in FP64. It produces exactly the same iteration counts and final residuals as cuSPARSE-, PETSc-, and Ginkgo-based FP64 CG, so the comparison isolates the precision controller from the execution engine. Fig. 7 shows that the learned controller improves time-to-solution while preserving FP64-level convergence quality. For readability, Fig. 7 visualizes 21 representative cases, but all aggregate values

Figure 7: Mixed-precision ablation on CG against EvoCGFP64, the double-precision solver on the identical execution substrate. (a) Speedup of EvoCG-Mixed over EvoCG-FP64. (b) Increase in iterations (Mixed vs. FP64). (c) Precision distribution per matrix. (d) Comparison of the final FP64-recomputed true-residual norms of EvoCG-Mixed and EvoCG-FP64 on a log10 scale. The dashed line denotes 𝑦 = 𝑥. below are computed over the full set of 47 matrices. Fig. 7 (a) reports an geometric-mean speedup of 1.28× over all 47 matrices, while showing 21 representative cases for readability; Fig. 7 (b) shows that this speedup is achieved with only an 11.4% average increase in iteration count; and Fig. 7 (c) shows substantial variation in precision mix across matrices, indicating matrix-dependent decisions rather than a fixed global policy, while more than 56%

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

Zhai et al.

of the SpMV computation is executed in low precision on average. For numerical validation, after each solve we recompute the true residual 𝑟 true = 𝑏 − 𝐴𝑥 in FP64 for both variants. Fig. 7(d) Mixed ∥ against log ∥𝑟 FP64 ∥ : the points closely track plots log10 ∥𝑟 true 2 10 true 2 the 𝑦 = 𝑥 line, with a geometric-mean Mixed/FP64 true-residual norm ratio of only 1.0108, confirming that EvoCG-Mixed reaches nearly the same final accuracy as EvoCG-FP64. The true-residual recomputation is performed outside the timed solver execution. Fig. 7(a)–(c) show that the benefit of dynamic precision control is dominated by a favorable group of 13 matrices, including jnlbrng1, minsurfo, and crystm02. On this set, mixed precision improves performance by 36.5% with only a 3.2% increase in iteration count, so the reduced SpMV cost translates almost directly into solverlevel speedup. A second group, including t3dl_e, bodyy4, ted_B, and Pres_Poisson, incurs a larger iteration increase of 31.8% on average but still delivers a 9% speedup, since the larger low-precision fraction outweighs the longer recurrence: iteration growth does not negate the benefit as long as per-iteration savings exceed the convergence penalty. Boundary cases are rare. thermal2 uses little low precision (86% FP64 and 14% FP32) yet still incurs a 26% iteration increase, leading to an 8% slowdown, while gridgena remains entirely in FP64 and thus exposes only the policy overhead, causing a modest 2% slowdown. Overall, dynamic precision control is robustly beneficial, and notably, no EvoCG-Mixed run diverges or misses the tolerance on any of the 47 matrices. In the worst cases, the policy autonomously retreats to a predominantly FP64 schedule, bounding the downside to the mild slowdowns above, roughly the policy overhead.

4.4

Analysis of Convergence Stability

A central question is whether EvoCG can exploit low precision while preserving CG’s essential convergence dynamics. Fig. 8 shows that EvoCG reaches the target tolerance and closely tracks the FP64 trajectory of the recursively updated residual. FP64: 15% FP32: 82% BF16: 3%

FP64: 18% FP32: 6% BF16: 76%

FP64: 66% FP32: 32% BF16: 2%

FP64: 4% FP32: 3% FP64: 5% BF16: 93% FP32: 2% BF16: 93%

FP64: 92% FP32: 5% BF16: 3%

FP64: 72% FP32: 15% BF16: 12%

FP64: 58% FP32: 42% BF16: 0% FP64: 60% FP32: 40% FP64: 43% BF16: 0% FP32: 57% BF16: 0%

FP64: 3% FP32: 0% BF16: 97%

FP64: 42% FP32: 58% BF16: 0%

FP64: 26% FP32: 19% BF16: 55%

FP64: 18% FP32: 20% BF16: 61%

Figure 8: Trajectories of the recursively updated relative residual for EvoCG-FP64 and EvoCG-Mixed on four representative matrices.

The first two cases, Dubcova2 and wathen100, show that EvoCGMixed can be nearly convergence-neutral. Dubcova2 converges in the same number of iterations as EvoCG-FP64, and wathen100 stays tightly aligned with only a moderate increase from 314 to 345 iterations even though BF16 dominates most later intervals. Both descend smoothly without late-stage stagnation or oscillation, so EvoCG can exploit low precision aggressively on tolerant problems at little or no penalty. The other two cases represent a more sensitive regime. On ted_B_unscaled, both solvers show the sharp non-monotone fluctuations typical of this problem, yet the mixed-precision run reaches the tolerance with only a small delay, 45 iterations versus 41. On bodyy4 the gap is larger, 268 versus 227, but the residual still decreases steadily. The penalty thus appears as reduced convergence efficiency rather than stagnation, breakdown, or loss of numerical control. The precision trajectories further show that the controller does not follow a fixed monotone schedule: wathen100, bodyy4 and ted_B_unscaled shift toward more BF16 in later stages, whereas Dubcova2 grows more conservative near the end. EvoCG thus adapts to the numerical tolerance of each solve rather than simply raising precision over time.

4.5

Cross-Scale Transfer of the Learned Precision Policy

Using 𝑝𝑘 as the policy state ties the learned policy to numerical distributions rather than to matrix-specific sparsity patterns, raising the possibility that a policy trained on one problem transfers to others. To test this, we train the controller only on mesh2e1, a tiny 306 × 306 matrix with 2018 nonzeros and under 15 minutes of training, then transfer it to all 47 unseen matrices without retraining. We do not claim universal transfer. The transferred policy outperforms EvoCG-FP64 on 23 of the 47 matrices. Across these 23 improved cases, it achieves a geometric-mean speedup of 1.21×. Among the remaining matrices, 15 converge without a performance gain and 9 fail to converge. Transfer is thus real across scales but bounded by problem family, since 𝑝𝑘 encodes numerical distributions shared within a family rather than across all matrices. We treat it as a zero-retraining fast path for problems near the training family, with matrix-specific training as the default and fallback whenever transfer fails. Fig. 9 reports five representative matrices spanning outcomes from better to worse than the matrix-specific policy: EvoCG-Mixed is the matrix-specific policy, EvoCG-Mixed (mesh2e1) the transferred policy, and EvoCG-FP64 the double-precision reference. On minsurfo, the transferred policy even improves over the matrixspecific one, cutting iterations from 103 to 97. On qa8fm the two are identical at 75 iterations, and on thermomech_dM they differ by only one iteration. On wathen120, transfer costs only a minor 354-to-360 increase. Dubcova1 is the only clear outlier, needing 176 iterations instead of 116.

4.6

Extension to Jacobi-Preconditioned CG

We further evaluate whether a policy trained on unpreconditioned CG remains effective when the solver dynamics are altered by preconditioning. We apply EvoCG to Jacobi-PCG on all 47 matrices

EvoCG : Self-Evolving Precision Control for Faster Conjugate Gradient

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

Figure 9: Representative convergent cases of cross-scale policy transfer from mesh2e1 to larger unseen matrices. without retraining, reusing the policy trained for the corresponding unpreconditioned system. The policy continues to control only the SpMV 𝑞𝑘 = 𝐴𝑝𝑘 , while the Jacobi application 𝑧𝑘 = 𝑀 −1𝑟𝑘 and all scalar operations remain in FP64. Both EvoCG-Mixed and EvoCGFP64 converge on all 47 matrices. EvoCG-Mixed achieves a 4.17× geometric-mean speedup over cuSPARSE-based Jacobi-PCG and a 1.15× speedup over EvoCG-FP64, with only a 4.99% average increase in iteration count. The latter comparison isolates the benefit of adaptive precision on the same preconditioned execution substrate, showing that the policies learned for CG transfer effectively to Jacobi-PCG without retraining.

4.7

Why Learned Precision Control?

We next ask whether EvoCtrl is necessary, or whether a simple rule-based scheduler can achieve comparable robustness and performance. To isolate the precision-control mechanism, we replace only EvoCtrl with a residual-threshold rule adapted from MilleMixed [46], while keeping the same Tile-CSC layout, mixed-precision SpMV kernel, and EvoRun runtime. Given the convergence tolerance 𝜖 = 10−10 , the rule assigns each tile a precision from the largest entry magnitude 𝑚𝑖 = max 𝑗 |𝑝𝑘(𝑖 ) [ 𝑗]| in its segment: BF16 if 𝑚𝑖 < 𝜖 × 10−1 , FP32 if 𝑚𝑖 ∈ [𝜖 × 10−1, 𝜖), and FP64 otherwise. Thus, the comparison attributes differences in convergence and execution time directly to the precision decisions.

EvoCG converges on all 47. Table 1 reports representative matrices covering both regimes. The conservative regime dominates the convergent cases shown in the table: on finan512, bcsstm39, cvxbqp1, jnlbrng1, bundle1, ted_B_unscaled, and cbuckle, the rule keeps most tiles in FP64, yielding iteration counts close to FP64 but leaving a substantial performance gap to EvoCG. In contrast, the aggressive regime appears on ted_B, where the rule requires 821 iterations compared with 493 for FP64, and more broadly in the 20 matrices where the rule fails to converge within the 10,000iteration limit. The table shows two representative non-convergent cases, Dubcova2 and Pres_Poisson. These results show that entry magnitude alone does not reliably capture tile-level numerical sensitivity. EvoCtrl optimizes time-to-solution rather than iteration count. Across the eight displayed matrices on which both methods converge, EvoCG is 1.36× faster than the rule on geometric-mean and is faster in every case. This remains true when EvoCG executes more iterations. On bcsstm39, for example, EvoCG requires 499 iterations instead of 415, yet reduces execution time from 18.37 ms to 12.76 ms. Overall, a scalar magnitude threshold cannot simultaneously expose substantial low-precision opportunities and preserve robust convergence across matrices. EvoCtrl instead learns statedependent, tile-wise precision decisions from the evolving Krylov process, allowing EvoCG to reduce per-iteration cost without entering the unstable regime of an overly aggressive fixed rule.

Table 1: Comparison of EvoCG, residual-threshold precision scheduling, and cuSPARSE FP64.

4.8 Iterations Matrix

EvoCG

Rule

Time (ms) FP64 EvoCG

Rule

FP64

finan512 58 (+ 1) 57 57 2.68 4.06 7.39 bcsstm39 499 (+ 84) 415 415 12.76 18.37 46.71 cvxbqp1 5095 (+ 14) 5083 (+ 2) 5081 161.05 281.30 583.76 jnlbrng1 156 (+ 2) 154 154 4.35 5.65 19.00 bundle1 264 (+ 49) 216 (+ 1) 215 11.50 15.98 27.11 ted_B 682 (+ 189) 821 (+ 328) 493 13.48 15.52 51.10 ted_B_unscaled 45 (+ 4) 41 41 1.13 1.28 5.05 cbuckle 4849 (+ 7) 4845 (+ 3) 4842 240.54 261.22 524.91 Dubcova2 193 10000 193 11.88 1571.47 25.51 Pres_Poisson 3428 (+ 993) 10000 2435 131.05 560.07 266.85

The rule exposes two undesirable operating regimes. It either retains mostly FP64 computation and provides limited speedup, or uses low precision aggressively and risks substantial iteration growth or convergence failure. Under the same 10,000-iteration limit, the rule converges on only 27 of the 47 matrices, whereas

Cost Analysis

While EvoCG improves solver performance through adaptive mixedprecision execution, it also incurs three extra costs: online policy execution, auxiliary Tile-CSC and precision metadata storage, and offline policy training. We analyze them separately to clarify the system overhead.

Policy cost. Fig. 10(a) reports the runtime overhead of embedding the learned policy into the CG loop. Across all 47 matrices, policy inference accounts for only 1.09% of end-to-end execution time on average, with a median of 0.60%. Moreover, 31 matrices remain below 1% overhead, and 38 remain below 2%. This result confirms that the lightweight two-layer MLP, fused device-side inference kernel, and profile-guided interval 𝐾 ★ together make neural precision control a practical in-loop mechanism rather than a solver-level bottleneck.

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

Memory cost. The memory overhead of Tile-CSC remains well controlled for the dominant portion of the workload set. As illustrated in Fig. 10 (b), on the 40 non-outlier matrices, Tile-CSC requires only 1.84× the memory of CSR overall. In this normal regime, the extra storage mainly comes from two deliberate choices: maintaining pre-quantized low-precision backups alongside the FP64 values, and using padding to regularize tile layout for warp-centric execution. This memory increase is necessary for the unified precision path of EvoCG in Sec. 3.3.2. By materializing multiple value streams offline, EvoCG can switch precision at runtime through direct value-stream selection, without incurring online quantization overhead or changing the traversal path. This design is important for preserving execution regularity and fully exposing the benefit of low-precision computation. A small number of outliers incur substantially higher memory overhead because their nonzeros are dispersed across many tiles, which magnifies tile-level structural and padding cost. Overall, Tile-CSC pays a controlled and necessary memory premium in exchange for efficient unified-kernel mixed-precision execution. Training cost. The training cost of EvoCG is entirely offline and amortized across future solves. In our experiments, the policy converges within 1000 episodes. Under this setting, training on a small matrix such as mesh2e1 takes less than 15 minutes. As discussed in Sec. 4.5, the mesh2e1-trained policy provides useful but nonuniversal cross-scale transfer, reducing the need for matrix-specific training on a subset of related problems. Matrix-specific training typically takes several hours and is therefore intended for workloads that repeatedly solve the same or closely related systems.

5

Related Work

Recent work accelerates CG and related Krylov solvers via tilegrained single-kernel GPU execution [46], adaptive mixed-precision PCG [22], RL-based precision tuning [12], and architectural support [37]. Broader mixed-precision Krylov methods explore inexact operators, residual correction, robust preconditioning, and nested solvers [26, 38, 45]. These studies establish the value of precision heterogeneity but treat it as an execution option, a safeguard, or an external tuning layer. EvoCG instead makes precision an online control variable derived from the solver state and consumed by tile-grained sparse execution in the loop. Table 2 summarizes this positioning. EvoCG also differs from the broader CG-acceleration literature. Prior work reduces communication, synchronization, and launch overheads through pipelining, communication-avoiding variants, and kernel fusion [2, 9, 10, 32, 36, 39]. Other work improves the sparse substrate through format redesign, vectorization, or load balancing [4, 14, 27, 28]. Mixed-precision numerical linear algebra

Figure 10: Cost analysis of EvoCG. (a) Policy inference time ratio across matrices. (b) Memory usage of Tile-CSC relative to CSR.

Zhai et al.

Table 2: Positioning of EvoCG relative to the closest mixedprecision CG systems.

Decision Mille-mixed [46] static rule Callipepla [37] fixed scheme RL tuning [12] RL (offline) Adaptive PCG [22] numerical rule EvoCG RL (in-loop)

Granularity

In-loop Closed Format loop co-design overhead

tile stream operator operator tile

✗ ✓ ✗ FPGA ISA ✗ ✗ partial ✗ ✓ ✓

✗ ✗ ✗ ✗ ✓

likewise trades arithmetic fidelity for performance under proper numerical control [1, 17, 25, 40]. EvoCG closes this loop, making precision selection a runtime solver action.

6

Conclusion

This paper presents EvoCG, a mixed-precision CG solver that treats tile-wise precision as a closed-loop control problem over the evolving Krylov state. By co-designing EvoCtrl, EvoSpMV, and EvoRun, EvoCG learns when reduced precision is numerically tolerable and translates it into end-to-end speedup. Across 47 SPD matrices from eight scientific domains, EvoCG delivers the best overall end-toend performance among all evaluated baselines and extends to Jacobi-PCG without retraining. More broadly, EvoCG suggests a path toward self-evolving HPC solvers that regulate their numerical behavior on modern low-precision hardware.

Acknowledgements This work was supported by the National Key Research and Development Program of China under Grant 2024YFB4505202, the National Natural Science Foundation of China under Grant 62632004, and the Hubei Provincial Natural Science Foundation of China under Grant 2026AFA002.

EvoCG : Self-Evolving Precision Control for Faster Conjugate Gradient

References [1] Ahmad Abdelfattah, Hartwig Anzt, Erik G Boman, Erin Carson, Terry Cojean, Jack Dongarra, Alyson Fox, Mark Gates, Nicholas J Higham, Xiaoye S Li, et al. 2021. A survey of numerical linear algebra methods utilizing mixed-precision arithmetic. The International Journal of High Performance Computing Applications 35, 4 (2021), 344–369. https://doi.org/10.1177/10943420211003313 [2] José I Aliaga, Joaquín Pérez, and Enrique S Quintana-Ortí. 2015. Systematic fusion of CUDA kernels for iterative sparse linear system solvers. In European Conference on Parallel Processing. Springer, 675–686. https://doi.org/10.1007/9783-662-48096-0_52 [3] Hartwig Anzt, Terry Cojean, Goran Flegar, Fritz Göbel, Thomas Grützmacher, Pratik Nayak, Tobias Ribizel, Yuhsiang Mike Tsai, and Enrique S Quintana-Ortí. 2022. Ginkgo: A modern linear operator algebra framework for high performance computing. ACM Transactions on Mathematical Software (TOMS) 48, 1 (2022), 1–33. https://doi.org/10.1145/3480935 [4] Hartwig Anzt, Terry Cojean, Chen Yen-Chen, Jack Dongarra, Goran Flegar, Pratik Nayak, Stanimire Tomov, Yuhsiang M Tsai, and Weichung Wang. 2020. Loadbalancing sparse matrix vector product kernels on GPUs. ACM Transactions on Parallel Computing (TOPC) 7, 1 (2020), 1–26. https://doi.org/10.1145/3380930 [5] Satish Balay, Shrirang Abhyankar, Mark F Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil Constantinescu, Lisandro Dalcin, Alp Dener, et al. 2022. PETSc: Portable, extensible toolkit for scientific computation. Astrophysics Source Code Library (2022), ascl–2210. [6] Nathan Bell and Michael Garland. 2009. Implementing sparse matrix-vector multiplication on throughput-oriented processors. In Proceedings of the conference on high performance computing networking, storage and analysis. 1–11. https: //doi.org/10.1145/1654059.1654078 [7] Amina Bouras and Valérie Frayssé. 2005. Inexact matrix-vector products in Krylov methods for solving linear systems: A relaxation strategy. SIAM J. Matrix Anal. Appl. 26, 3 (2005), 660–678. https://doi.org/10.1137/S0895479801384743 [8] Erin Carson and Nicholas J Higham. 2018. Accelerating the solution of linear systems by iterative refinement in three precisions. SIAM Journal on Scientific Computing 40, 2 (2018), A817–A847. [9] Erin C Carson. 2018. The adaptive 𝑠 -step conjugate gradient method. SIAM J. Matrix Anal. Appl. 39, 3 (2018), 1318–1338. https://doi.org/10.1137/16M1107942 [10] Erin C Carson, Miroslav Rozloznik, Zdenek Strakos, Petr Tichy, and Miroslav Tůma. 2018. The numerical stability analysis of pipelined conjugate gradient methods: Historical context and methodology. SIAM Journal on Scientific Computing 40, 5 (2018), A3549–A3580. https://doi.org/10.1137/16M1103361 [11] Ru-Shan Chen, Edward Kai-Ning Yung, Chi H Chan, Dao Xiang Wang, and Da Gang Fang. 2002. Application of the SSOR preconditioned CG algorithm to the vector FEM for 3D full-wave analysis of electromagnetic-field boundary-value problems. IEEE Transactions on Microwave Theory and Techniques 50, 4 (2002), 1165–1172. http://doi.org/10.1109/22.993420 [12] Xinye Chen. 2025. Mixed-precision conjugate gradient solvers with RL-driven precision tuning. arXiv preprint arXiv:2504.14268 (2025). https://doi.org/10.48550/ arXiv.2504.14268 [13] Chin-Hsiang Cheng, Yu-Xian Huang, Shun-Chih King, Chun-I Lee, and ChihHsing Leu. 2014. CFD (computational fluid dynamics)-based optimal design of a micro-reformer by integrating computational fluid dynamics code using a simplified conjugate-gradient method. Energy 70 (2014), 355–365. https: //doi.org/10.1016/j.energy.2014.04.005 [14] Kazem Cheshmi, Zachary Cetinic, and Maryam Mehri Dehnavi. 2022. Vectorizing sparse matrix computations with partially-strided codelets. In SC22: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 1–15. https://doi.org/10.1109/SC41404.2022.00037 [15] Jack Choquette. 2023. NVIDIA Hopper H100 GPU: Scaling performance. IEEE Micro 43, 3 (2023), 9–17. https://doi.org/10.1109/MM.2023.3256796 [16] Jack Choquette, Wishwesh Gandhi, Olivier Giroux, Nick Stam, and Ronny Krashinsky. 2021. NVIDIA A100 tensor core GPU: Performance and innovation. IEEE Micro 41, 2 (2021), 29–35. https://doi.org/10.1109/MM.2021.3061394 [17] Michael A Clark, Ronald Babich, Kipton Barros, Richard C Brower, and Claudio Rebbi. 2010. Solving lattice QCD systems of equations using mixed precision solvers on GPUs. Computer Physics Communications 181, 9 (2010), 1517–1528. https://doi.org/10.1016/j.cpc.2010.05.002 [18] Timothy A Davis and Yifan Hu. 2011. The University of Florida sparse matrix collection. ACM Transactions on Mathematical Software (TOMS) 38, 1 (2011), 1–25. https://doi.org/10.1145/2049662.2049663 [19] Jack Dongarra, John Gunnels, Harun Bayraktar, Azzam Haidar, and Dan Ernst. 2024. Hardware trends impacting floating-point computations in scientific applications. arXiv preprint arXiv:2411.12090 (2024). https://doi.org/10.48550/arXiv. 2411.12090 [20] Gene Howard Golub, Tony F Chan, and Raymond H Chan. 1997. Iterative Methods in Scientific Computing. Springer. [21] Anne Greenbaum. 1989. Behavior of slightly perturbed Lanczos and conjugategradient recurrences. Linear Algebra Appl. 113 (1989), 7–63. https://doi.org/10. 1016/0024-3795(89)90285-1

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

[22] Yichen Guo, Eric de Sturler, and Tim Warburton. 2025. An adaptive mixed precision and dynamically scaled preconditioned conjugate gradient algorithm. arXiv preprint arXiv:2505.04155 (2025). https://doi.org/10.48550/arXiv.2505.04155 [23] Suyog Gupta, Ankur Agrawal, Kailash Gopalakrishnan, and Pritish Narayanan. 2015. Deep learning with limited numerical precision. In International conference on machine learning. PMLR, 1737–1746. [24] Magnus R Hestenes, Eduard Stiefel, et al. 1952. Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards 49, 6 (1952), 409–436. https://doi.org/10.6028/JRES.049.044 [25] Nicholas J Higham and Theo Mary. 2022. Mixed precision algorithms in numerical linear algebra. Acta Numerica 31 (2022), 347–414. https://doi.org/10.1017/ S0962492922000022 [26] Takuya Ina, Yasuhiro Idomura, Toshiyuki Imamura, and Naoyuki Onodera. 2023. A new data conversion method for mixed precision Krylov solvers with FP16/BF16 Jacobi preconditioners. In Proceedings of the International Conference on High Performance Computing in Asia-Pacific Region. 29–34. https://doi.org/10.1145/ 3578178.3578222 [27] Weifeng Liu and Brian Vinter. 2015. CSR5: An efficient storage format for crossplatform sparse matrix-vector multiplication. In Proceedings of the 29th ACM on International Conference on Supercomputing. 339–350. https://doi.org/10.1145/ 2751205.2751209 [28] Duane Merrill and Michael Garland. 2016. Merge-based parallel sparse matrixvector multiplication. In SC’16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 678–689. https://doi.org/10.1109/SC.2016.57 [29] Paulius Micikevicius, Sharan Narang, Jonah Alben, Gregory Diamos, Erich Elsen, David Garcia, Boris Ginsburg, Michael Houston, Oleksii Kuchaiev, Ganesh Venkatesh, et al. 2018. Mixed precision training. In International Conference on Learning Representations. https://doi.org/10.48550/arXiv.1710.03740 [30] Richard Tran Mills, Mark F Adams, Satish Balay, Jed Brown, Alp Dener, Matthew Knepley, Scott E Kruger, Hannah Morgan, Todd Munson, Karl Rupp, et al. 2021. Toward performance-portable PETSc for GPU-based exascale systems. Parallel Comput. 108 (2021), 102831. https://doi.org/10.1016/j.parco.2021.102831 [31] Maxim Naumov, L Chien, Philippe Vandermersch, and Ujval Kapasi. 2010. cuSPARSE library. In GPU Technology Conference, Vol. 12. [32] Karl Rupp, Josef Weinbub, Ansgar Jüngel, and Tibor Grasser. 2016. Pipelined iterative solvers with kernel fusion for graphics processing units. ACM Transactions on Mathematical Software (TOMS) 43, 2 (2016), 1–27. https://doi.org/10.1145/2907944 [33] Yousef Saad. 2003. Iterative Methods for Sparse Linear Systems. SIAM. [34] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. 2017. Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347 (2017). https://doi.org/10.48550/arXiv.1707.06347 [35] Valeria Simoncini and Daniel B Szyld. 2003. Theory of inexact Krylov subspace methods and applications to scientific computing. SIAM Journal on Scientific Computing 25, 2 (2003), 454–477. https://doi.org/10.1137/S1064827502406415 [36] Edgar Solomonik, Erin Carson, Nicholas Knight, and James Demmel. 2017. Tradeoffs between synchronization, communication, and computation in parallel linear algebra computations. ACM Transactions on Parallel Computing (TOPC) 3, 1 (2017), 1–47. https://doi.org/10.1145/2897188 [37] Linghao Song, Licheng Guo, Suhail Basalama, Yuze Chi, Robert F Lucas, and Jason Cong. 2023. Callipepla: Stream centric instruction set and mixed precision for accelerating conjugate gradient solver. In Proceedings of the 2023 ACM/SIGDA International Symposium on Field Programmable Gate Arrays. 247–258. https: //doi.org/10.1145/3543622.3573182 [38] Kengo Suzuki and Takeshi Iwashita. 2025. A nested Krylov method using halfprecision arithmetic. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. 711–727. https: //doi.org/10.1145/3712285.3759807 [39] James D Trotter, Sinan Ekmekçibaşı, Doğan Sağbili, Johannes Langguth, Xing Cai, and Didem Unat. 2025. CPU- and GPU-initiated communication strategies for conjugate gradient methods on large GPU clusters. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. 298–315. https://doi.org/10.1145/3712285.3759774 [40] Jasper Van Den Eshof and Gerard LG Sleijpen. 2004. Inexact Krylov subspace methods for linear systems. SIAM J. Matrix Anal. Appl. 26, 1 (2004), 125–153. https://doi.org/10.1137/S0895479802403459 [41] Jose O Vargas and Ricardo Adriano. 2022. Subspace-based conjugate-gradient method for solving inverse scattering problems. IEEE Transactions on Antennas and Propagation 70, 12 (2022), 12139–12146. https://doi.org/10.1109/TAP.2022. 3211738 [42] Bastien Vieublé. 2022. Mixed precision iterative refinement for the solution of large sparse linear systems. Ph. D. Dissertation. INP Toulouse. [43] Samuel Williams, Andrew Waterman, and David Patterson. 2009. Roofline: An insightful visual performance model for multicore architectures. Commun. ACM 52, 4 (2009), 65–76. https://doi.org/10.1145/1498765.1498785 [44] Zhengpeng Xie, Qiang Zhang, Fan Yang, Marco Hutter, and Renjing Xu. 2024. Simple policy optimization. arXiv preprint arXiv:2401.16025 (2024). https://doi. org/10.48550/arXiv.2401.16025

PPoPP’27, March 20–24, 2027, Salt Lake City, Utah

[45] Ichitaro Yamazaki, Erin Carson, and Brian Kelley. 2022. Mixed precision 𝑠 -step conjugate gradient with residual replacement on GPUs. In 2022 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 886–896. https: //doi.org/10.1109/IPDPS53621.2022.00091

Zhai et al.

[46] Dechuang Yang, Yuxuan Zhao, Yiduo Niu, Weile Jia, En Shao, Weifeng Liu, Guangming Tan, and Zhou Jin. 2024. Mille-feuille: A tile-grained mixed precision single-kernel conjugate gradient solver on GPUs. In SC24: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 1–16. https://doi.org/10.1109/SC41406.2024.00064

Related documents

• MANOJAVAM: A Scalable, Unified FPGA Accelerator for Matrix Multiplication and Singular Value Decomposition in Principal Component Analysis#155256
• Nora: Normalized Orthogonal Row Alignment for Scalable Matrix Optimizer#155302
• Mixed-Precision Computing for Scientific Discovery: Formats, Co-Design, and Responsible Approximation#1132338
• Exposing Embedded Matrix Multiplication with polyhedral transformations for CGRA compilation#1400900
• Hierarchical Muon: Tiled Newton-Schulz Updates for Efficient Muon Optimization#310817