{"id":"5cb0fb0a-cb54-4ecb-8694-8210b1d86104","arxiv_id":"1908.02964","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":4.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"A straightforward extension of Bayesian conjugate gradient solves several related linear systems simultaneously by treating them as one block-diagonal system with a joint prior.","lead":"This short discussion piece proposes a way to solve multiple related sets of linear equations together in one Bayesian framework, so that solving one can inform the others. It is a comment on a companion paper about a Bayesian version of the conjugate gradient method.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Claimed reduction to independent BayesCG with block-diagonal prior is false; cross-system dependence arises from joint observations, not solely from Σ0.","rationale":"After carefully checking the mathematics, the reader's identified weakest assumption—the search-direction recurrence—does not actually pose a correctness risk. The recurrence is the standard CG recurrence in the M-inner product on the residual space; with M = AΣ0A^T it generates M-orthogonal search directions. The original BayesCG paper's proof carries over directly to the larger block system. However, a more serious and concrete error is present: the claimed reduction to J independent pairs of BayesCG for block-diagonal Σ0 is false. The joint algorithm's observations are scalar combinations of all systems' right-hand sides, and the recurrence coefficient is a sum over systems, so systems are coupled even when the prior covariance has no off-diagonal terms. This is demonstrated by the simple example above. This error undermines the paper's characterization of the source of cross-system dependence (Remark 1) and its 'strict generalisation' statement (Remark 3), but it does not invalidate the core simultaneous-solving algorithm. The paper should be revised to correct these remarks and clarify that coupling is inherent to the joint observation model; the prior Σ0 modulates rather than creates the coupling. Therefore the existing CONDITIONAL verdict remains appropriate, altered only by the specific reason.","tokens_in":4057,"tokens_out":37511,"duration_ms":324805,"concrete_test":"Run the multi-system BayesCG algorithm on J=2, d=2 with A1=diag(2,1), A2=diag(1,2), x0=0, Σ0=I_4, and b1=(1,2), b2=(3,4). Record the first iterate for system 1 and compare to the first iterate of BayesCG applied to system 1 alone. Then rerun with b2=(4,4). If the joint first iterate differs from the independent iterate, and if changing b2 changes the joint estimate for system 1, then Remark 3's reduction claim is false and cross-system dependence occurs independently of Σ0.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The reader's concern about the unproved search-direction recurrence is not the weakest point. The recurrence \\tilde{s}_m = r_{m-1} - ⟨s_{m-1}, r_{m-1}⟩_M s_{m-1} with M = AΣ0A^T is exactly the conjugate-gradient recurrence for M z = r0, so the s_m are M-orthogonal without further assumptions. The real problem is Remarks 1 and 3. When Σ0 is block-diagonal, the algorithm does not reduce to J separate BayesCG instances: the scalar observation y_m = s_m^T b = Σ_j s_{m,j}^T b_j and the recurrence coefficient ⟨s_{m-1}, r_{m-1}⟩_M = Σ_j s_{m-1,j}^T M_j r_{m-1,j} are sums over systems, so search directions and posterior updates couple systems even with an independent prior. Hence the estimator for x_j depends on b_{j'} for j' ≠ j even when Σ0 is block-diagonal, contradicting Remark 1's claim that dependence comes from Σ0. For example, with A1=diag(2,1), A2=diag(1,2), b1=(1,2), b2=(3,4), x0=0, Σ0=I_4, the first joint iterate for system 1 is (20/27,20/27), whereas independent BayesCG gives (5/4,5/4); changing b2 to (4,4) changes the joint iterate to (37/44,37/44). Thus the block-diagonal case is not a strict reduction, and the proposed transfer mechanism is confounded by algorithmic coupling.","agreement_with_reader":"disagree"},"referee_report":{"model":"deepseek-v4-flash","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.","tokens_in":4332,"tokens_out":13378,"duration_ms":125164,"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":[{"comment":"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.","section":"A generalisation to multiple linear systems, Remark 3"},{"comment":"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.","section":"A generalisation to multiple linear systems, Remark 1"}],"minor_comments":[{"comment":"There is a typo in 'p aper' in both the abstract and the opening sentence; it should be 'paper'.","section":"Abstract and opening paragraph"},{"comment":"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.","section":"A generalisation to multiple linear systems, Remark 2"},{"comment":"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.","section":"A generalisation to multiple linear systems, search-direction definition"},{"comment":"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.","section":"A generalisation to multiple linear systems, notation"}],"recommendation":"major_revision","confidential_remarks":null},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The paper's actual contribution is a proposed extension of BayesCG to solve several related linear systems at once by stacking them into a block system with a joint prior, plus a sensible discussion of prior elicitation and separable covariances as a transfer-learning mechanism. It is clearly written and honestly labels its own limitations, for example the lack of a sequential variant. The search-direction recurrence the reader worried about is actually the standard CG recurrence for M-orthogonality with M = AΣ0A^T; that part is fine.\n\nThe real problem is Remarks 1 and 3. When Σ0 is block-diagonal, the algorithm does not reduce to J independent BayesCG instances. The scalar observation y_m = s_m^T b and the recurrence coefficient ⟨s_{m-1}, r_{m-1}⟩_M are sums over the J systems, so search directions and posterior updates couple the systems even with an independent prior. A concrete check: take A1=diag(2,1), A2=diag(1,2), b1=(1,2), b2=(3,4), x0=0, Σ0=I_4. The first joint iterate for system 1 is (20/27,20/27), whereas independent BayesCG gives (5/4,5/4). Change b2 to (4,4) and the joint iterate becomes (37/44,37/44). So the estimator for x1 depends on b2 even when the prior is block-diagonal. Remark 3's strict-generalisation claim is false, and Remark 1's attribution of all cross-system dependence to Σ0 is wrong. The coupling is real but it comes from the algorithmic observation model, not solely from the prior.\n\nThis does not kill the transfer-learning idea, but it means the paper's stated mechanism is incorrect. The authors make no numerical experiments, so the claimed accuracy and uncertainty benefits are untested. The derivation of the posterior updates is standard and sound; the flaw is in the interpretation of the block-diagonal special case.\n\nFor a contributed discussion this is a reasonable prompt, but the central technical claim overreaches. The paper deserves serious peer review because referees should catch this and push for a corrected presentation; the idea itself is worth exploring, but only after the coupling is properly understood and demonstrated.\n\nFor you: worth a quick read if you work on probabilistic linear solvers, but do not cite it as a strict generalisation until the reduction claim is fixed.","headline":"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.","tokens_in":4902,"tokens_out":1923,"would_cite":false,"duration_ms":21697,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62F15","65F10"],"pacs":[],"model":"deepseek-v4-flash","headline":"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.","keywords":["Bayesian conjugate gradient","multi-system linear solver","probabilistic numerics","prior covariance","Kronecker product","transfer learning","Krylov subspace recycling"],"falsifier":"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.","tokens_in":3740,"feed_emoji":"🔗","tokens_out":8860,"duration_ms":90653,"temperature":0.7,"pith_summary":"This discussion of the Bayesian conjugate gradient method makes two points. First, prior choice for linear solvers should exploit known structure of the matrix, such as sparsity patterns inherited from a finite-element mesh, rather than only placing a prior on the solution vector. Second, BayesCG extends naturally to multiple related systems: stack the systems into one block-diagonal system and put a joint prior on the combined solution. In this multi-system BayesCG, the estimator for one system depends on the other systems through the prior covariance, which can improve accuracy and uncertainty quantification. The authors present the algorithm and special prior forms, but include no numerical experiments.","feed_headline":"Joint prior lets BayesCG solve related linear systems together","feed_subtitle":"Stack several systems into one block problem; each solution draws on the others through the prior covariance.","key_machinery":"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$.","core_discovery":"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.","pith_inferences":["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."],"forward_implications":["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."],"supporting_citations":[{"why":"Supplies the single-system Bayesian conjugate gradient algorithm whose equations and orthogonality requirement are extended to the block-diagonal system.","marker":"[2]"},{"why":"Justifies putting the prior on the solution rather than on the inverse of the matrix, the modelling choice the joint prior inherits.","marker":"[6]"},{"why":"Documents the need for solving several related linear systems and provides the non-Bayesian Krylov recycling alternative the authors build on.","marker":"[5]"},{"why":"Introduces the separable covariance construction $\\mathbf{B}\\otimes\\boldsymbol{\\Sigma}_0$ used as the interpretable prior for transfer across systems.","marker":"[8]"},{"why":"Motivates the broader goal of propagating uncertainty through chains of coupled computations, which the multi-system extension serves.","marker":"[7]"}],"fun_headline_variants":["Joint prior lets BayesCG solve multiple linear systems at once","BayesCG extended: one joint prior, many related systems","Solving related linear systems together with BayesCG prior","Block-diagonal prior generalizes BayesCG to multiple systems","BayesCG variant uses joint prior to solve systems in parallel"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["Joint prior lets BayesCG solve multiple linear systems at once","BayesCG extended: one joint prior, many related systems","Solving related linear systems together with BayesCG prior","Block-diagonal prior generalizes BayesCG to multiple systems","BayesCG variant uses joint prior to solve systems in parallel"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000194,"raw_usage":{"total_tokens":1316,"prompt_tokens":869,"completion_tokens":447,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":485,"completion_tokens_details":{"reasoning_tokens":366}},"tokens_in":485,"tokens_out":447,"duration_ms":5111,"temperature":1.0,"reasoning_tokens":366,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T14:28:17.845820+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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.","supporting_citations":[{"cited_title":"Cockayne, C","cited_arxiv_id":null,"evidence_quote":"Supplies the single-system Bayesian conjugate gradient algorithm whose equations and orthogonality requirement are extended to the block-diagonal system."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Justifies putting the prior on the solution rather than on the inverse of the matrix, the modelling choice the joint prior inherits."},{"cited_title":"Krylov Subspace Recycling for Fast Iterative Least-Squares in Machine Learning","cited_arxiv_id":"1706.00241","evidence_quote":"Documents the need for solving several related linear systems and provides the non-Bayesian Krylov recycling alternative the authors build on."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces the separable covariance construction $\\mathbf{B}\\otimes\\boldsymbol{\\Sigma}_0$ used as the interpretable prior for transfer across systems."},{"cited_title":"Hennig, M","cited_arxiv_id":null,"evidence_quote":"Motivates the broader goal of propagating uncertainty through chains of coupled computations, which the multi-system extension serves."}],"review_version":1}