REVIEW 2 major objections 4 minor 8 references
Contributed Discussion of "A Bayesian Conjugate Gradient Method"
T0 review · 2 major / 4 minor · reviewed 2026-08-14 · deepseek-v4-flash
Pith's one-line read The BayesCG algorithm can be extended to solve several related linear systems at once by stacking them into a block-diagonal system and coupling the solutions through a joint Gaussian prior.
desk verdict Proposed multi-system BayesCG is a natural idea, but the paper's own reduction claim is false and the transfer mechanism is confounded by algorithmic coupling. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The central object is the block-diagonal coefficient matrix $\mathbf{A}=\operatorname{BlockDiag}[\mathbf{A}_1,\dots,\mathbf{A}_J]$ together with the joint prior $\mathcal{N}(\mathbf{x}_0,\boldsymbol{\Sigma}_0)$, whose cross-system off-diagonal blocks carry the information transfer. The algorithm's operative mechanism is the requirement that search directions be orthogonal in the $\mathbf{A}\boldsymbol{\Sigma}_0\mathbf{A}^\top$ inner product; this is what lets a direction probed in one system inform the estimate in another. The separable prior $\boldsymbol{\Sigma}_0=\mathbf{B}\otimes\boldsymbol{\Sigma}_0$ is presented as the practical way to encode task correlation through $\mathbf{B}$ and per-system covariance through $\boldsymbol{\Sigma}_0$.
What would settle it
Take a small block system, say two $2\times2$ positive-definite systems, run the proposed recurrence for a few iterations with a generic positive-definite $\boldsymbol{\Sigma}_0$, and compute the $\mathbf{A}\boldsymbol{\Sigma}_0\mathbf{A}^\top$-inner products between distinct search directions; any nonzero off-diagonal value would show the directions are not orthogonal and would falsify the claim that this is a conjugate gradient method.
Extended reading notes
Core claim
The paper's central claim is that a joint Gaussian prior over the concatenated solutions of $J$ linear systems, combined with the same projection-by-projection observation model as in BayesCG, yields a genuine conjugate-gradient algorithm on the block-diagonal system $\mathbf{A}=\operatorname{BlockDiag}[\mathbf{A}_1,\dots,\mathbf{A}_J]$. Conditioning on the observed projections gives a Gaussian posterior with the same algebraic form as the single-system case but with larger vectors and matrices. The search directions are required to be $\mathbf{A}\boldsymbol{\Sigma}_0\mathbf{A}^\top$-orthogonal, and the paper asserts that the recurrence $\tilde{\mathbf{s}}_m=\mathbf{r}_{m-1}-\langle \mathbf{s}_{m-1},\mathbf{r}_{m-1}\rangle_{\mathbf{A}\boldsymbol{\Sigma}_0\mathbf{A}^\top}\mathbf{s}_{m-1}$ produces them. From this, the estimator for any one system is influenced by the other systems' matrices and right-hand sides through $\boldsymbol{\Sigma}_0$. The paper also notes that a block-diagonal $\boldsymbol{\Sigma}_0$ recovers $J$ independent BayesCG runs, so the multi-system version is a strict generalisation.
Load-bearing premise
The whole construction stands on the unproved assertion that the recurrence for the search directions produces directions that are mutually perpendicular in the geometry induced by the block matrix and the prior covariance; if that orthogonality fails, the method is not a conjugate gradient algorithm and the cross-system dependence is not achieved.
Editorial extensions
If this is right
- The multi-system estimator for $\mathbf{x}^*_j$ depends on data from the other systems whenever $\boldsymbol{\Sigma}_0$ has off-block-diagonal mass, so prior specification controls the strength of transfer.
- A block-diagonal $\boldsymbol{\Sigma}_0$ reproduces $J$ independent BayesCG runs, making the multi-system algorithm a strict generalisation rather than a different method.
- The separable form $\boldsymbol{\Sigma}_0=\mathbf{B}\otimes\boldsymbol{\Sigma}_0$ gives an interpretable separation between cross-task correlation and within-task covariance, so multi-output Gaussian process priors can be reused.
- Because one projection is observed from every system at each iteration, the algorithm solves the systems simultaneously; a sequential variant would require additional methodology.
- The cost rises to matrix-vector products with $dJ\times dJ$ matrices, so the method is attractive only when cross-system sharing buys enough accuracy or calibration.
Reading between the lines
- If the orthogonality recurrence holds, the method is effectively a Bayesian transfer-learning solver; a natural test is to compare its accuracy and uncertainty calibration against independently solved systems on a family of nearby linear systems from a Gaussian-process model-selection problem.
- Nothing in the construction ties $\boldsymbol{\Sigma}_0$ to the block structure, so any positive-definite cross-covariance from multi-output Gaussian processes or vector-valued kernels could be plugged in; the paper only illustrates the separable case.
- The sequential limitation suggests a natural extension: use the posterior from solving $\mathbf{A}_1$ as the prior for $\mathbf{A}_2$, which would turn the joint construction into an online transfer-learning method the paper explicitly leaves open.
- The prior-elicitation discussion implies that encoding spectral or sparsity information in $\boldsymbol{\Sigma}_0$ might accelerate convergence, but the paper does not demonstrate this; a synthetic test with matrices of known spectrum would settle it.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. This contributed discussion paper first comments on prior specification for Bayesian numerical linear solvers, arguing that prior knowledge about the matrix A (e.g., sparsity, spectrum, or conditioning) is natural and likely beneficial. The main contribution is a proposed generalization of the BayesCG algorithm of Cockayne et al. to solve J related linear systems Aj x*_j = b_j simultaneously. The authors form a block-diagonal operator A = BlockDiag[A1,...,AJ], place a joint Gaussian prior on the concatenated solution vector, and apply the BayesCG conditioning formulas to the larger system. They give the posterior mean and covariance updates, define search directions via the recurrence \tilde{s}_m = r_{m-1} - <s_{m-1}, r_{m-1}>_{A Sigma0 A^T} s_{m-1}, and make four remarks about cross-system dependence, separable priors, a claimed reduction to J independent BayesCG runs when Sigma0 is block-diagonal, and the simultaneous nature of the resulting algorithm.
Significance. If correct, the proposed multi-system extension would be a useful framework for transfer learning in probabilistic numerical linear algebra, aligning with current interest in reusing information across related computational tasks. The paper also raises a legitimate and underexplored question about how to elicit priors that reflect structural knowledge of the operator A. The algebraic conditioning steps are standard and appear correct, and the search-direction recurrence is a recognizable conjugate-gradient construction. The main novelty, however, is the claim that cross-system dependence can be controlled through Sigma0 and that a block-diagonal prior yields J independent BayesCG instances; this claim is demonstrably incorrect as stated, which undercuts the paper's central interpretation.
major comments (2)
- [A generalisation to multiple linear systems, Remark 3] The statement that a block-diagonal prior Sigma0 = BlockDiag[Sigma0,1,...,Sigma0,J] reduces the multi-system method to J separate instances of BayesCG is false. The scalar observation y_m = s_m^T b = sum_j s_{m,j}^T b_j and the recurrence coefficient <s_{m-1}, r_{m-1}>_M = sum_j s_{m-1,j}^T M_j r_{m-1,j}, with M_j = A_j Sigma0_j A_j^T, are sums over all systems, so the search directions and posterior updates couple the systems even when Sigma0 has no off-diagonal blocks. For example, take A1=diag(2,1), A2=diag(1,2), b1=(1,2), x0=0, and Sigma0=I4. The first joint iterate for system 1 is (20/27,20/27) when b2=(3,4), but changing b2 to (4,4) changes the first joint iterate to (37/44,37/44); independent BayesCG would give (5/4,5/4) in both cases. Thus the block-diagonal case is not a strict reduction, and the claimed reduction to J independent instances is incorrect.
- [A generalisation to multiple linear systems, Remark 1] The claim that cross-system dependence 'will come from the matrix Sigma0' is not supported. In the counterexample above with Sigma0 = I4, the estimator for x1 depends on b2 even though the prior is block-diagonal with independent blocks. The dependence arises from the global projection observations y_m and from the coupled search-direction recurrence; Sigma0 is not the only source of coupling. The remark should be revised to acknowledge that the algorithmic construction itself induces dependence across systems, regardless of the prior.
minor comments (4)
- [Abstract and opening paragraph] There is a typo in 'p aper' in both the abstract and the opening sentence; it should be 'paper'.
- [A generalisation to multiple linear systems, Remark 2] The notation 'Sigma0 = B ⊗ Sigma0' reuses Sigma0 for both the joint prior covariance and the per-system covariance, which is confusing; please use distinct symbols, e.g., Sigma0 = B ⊗ Sigma_tilde0.
- [A generalisation to multiple linear systems, search-direction definition] The recurrence \tilde{s}_m = r_{m-1} - <s_{m-1}, r_{m-1}>_{A Sigma0 A^T} s_{m-1} is stated without proof or reference. It is a standard conjugate-gradient construction, but for a self-contained discussion the authors should either provide a short justification or explicitly refer to the corresponding derivation in the original BayesCG paper.
- [A generalisation to multiple linear systems, notation] The line 'A x* = b where x = ((x*_1)^T, ..., (x*_J)^T)^T' appears to have a typo (extra asterisk or missing subscript); please clarify the definition of the concatenated vector.
Circularity Check
No significant circularity; the multi-system BayesCG is a direct application of BayesCG to a block-diagonal system, with no fitted inputs and no load-bearing self-citations.
full rationale
The paper's extension consists in embedding J related linear systems into a single block-diagonal system A = BlockDiag[A1,...,AJ] and then applying the BayesCG update equations verbatim to the larger vectors and matrices. This is a direct application of an existing algorithm to a new system, not a derivation that assumes its own conclusion. The prior covariance Sigma0 is an input to the algorithm, and the claimed cross-system dependence is a mathematical consequence of the update equations; no parameter is fitted to a subset of data and then renamed a prediction. The only self-citations, such as [8] for separable covariance functions, are motivational and are not load-bearing for the algorithm's derivation. The unproved recurrence for the search directions is presented as a straightforward analogue of BayesCG; even if one doubted its sufficiency, that would be a correctness or completeness issue, not circularity. Likewise, Remark 3's claim that a block-diagonal Sigma0 reduces the method to J independent BayesCG instances may be incorrect, because the joint normalization and summed observations can couple the systems even with an independent prior; however, an incorrect reduction claim does not make the conclusion equivalent to its inputs by construction. Therefore the derivation chain is self-contained and not circular.
Assumptions & free parameters
assumptions (3)
- domain assumption The posterior distribution for a linear system with Gaussian prior and noise-free projections is Gaussian with mean and covariance as given in Cockayne et al. [2].
- ad hoc to paper The search directions defined by the recurrence \tilde{s}_m = r_{m-1} - <s_{m-1}, r_{m-1}>_{A Sigma0 A^T} s_{m-1} are A Sigma0 A^T-orthogonal for the block system.
- domain assumption The matrix A is such that A Sigma0 A^T defines a valid inner product (e.g., A is symmetric positive definite and Sigma0 is positive semidefinite).
Cite this review
Pith. "Pith review of Contributed Discussion of "A Bayesian Conjugate Gradient Method"." pith.science (2026). https://pith.science/paper/4WRPVU3U
@misc{pith2026190802964,
author = {Pith},
title = {Pith review of: Contributed Discussion of "A Bayesian Conjugate Gradient Method"},
year = {2026},
howpublished = {\url{https://pith.science/paper/4WRPVU3U}},
note = {Machine review of arXiv:1908.02964}
}
read the original abstract
We would like to congratulate the authors of "A Bayesian Conjugate Gradient Method" on their insightful paper, and welcome this publication which we firmly believe will become a fundamental contribution to the growing field of probabilistic numerical methods and in particular the sub-field of Bayesian numerical methods. In this short piece, which will be published as a comment alongside the main paper, we first initiate a discussion on the choice of priors for solving linear systems, then propose an extension of the Bayesian conjugate gradient (BayesCG) algorithm for solving several related linear systems simultaneously.
Figures
Reference graph
Works this paper leans on
- [1]
-
[2]
J. Cockayne, C. J. Oates, I. C. F. Ipsen, and M. Girolami. A Bayesian conjugate gradient method (with discussion). Bayesian Analysis , 2019
work page 2019
-
[3]
Probabilistic Numerical Methods for Partial Differential Equations and Bayesian Inverse Problems
J. Cockayne, C. J. Oates, T. Sullivan, and M. Girolami. Pr obabilistic meshless methods for partial differential equations and Bayesian inverse proble ms. arXiv:1605.07811, 2016
work page Pith review arXiv 2016
-
[4]
T. A. Davis and Y. Hu. The University of Florida sparse mat rix collection. ACM Transactions on Mathematical Software , 38(1), 2011
work page 2011
-
[5]
Krylov Subspace Recycling for Fast Iterative Least-Squares in Machine Learning
F. de Roos and P. Hennig. Krylov subspace recycling for fa st iterative least-squares in machine learning. arXiv:1706.00241, 2017
work page Pith review arXiv 2017
-
[6]
P. Hennig. Probabilistic interpretation of linear solv ers. SIAM Journal on Optimization , 25, 2015
work page 2015
- [7]
-
[8]
X. Xi, F-X. Briol, and M. Girolami. Bayesian quadrature f or multiple related integrals. In International Conference on Machine Learning, PMLR 80 , pages 5369–5378, 2018. 4
work page 2018
Reviewed August 14, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.