{"id":"8c231e0b-5d4a-44f2-801b-961468552f78","arxiv_id":"2501.03377","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":3.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A pseudoinverse preconditioner built from the Kronecker structure of the finite-difference Laplacian makes PCG solve cuboid Poisson problems in two or three iterations.","lead":"This paper tests a fast math trick called a pseudoinverse preconditioner for solving the Poisson equation on rectangular grids used in ferroelectric material simulations. The trick makes the conjugate gradient solver finish in just a few steps, and the paper shows how to build and apply it efficiently.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The few-iteration convergence claim lacks a defined stopping criterion and reported iteration counts, so 'two or three sufficient' is not reproducible from the data.","rationale":"After reading the full text, I find the mathematical construction of the pseudoinverse preconditioner to be sound: for a Kronecker-sum Laplacian with orthonormal 1D eigenbases, Algorithms 6 and 7 apply the exact Moore-Penrose pseudoinverse, so PCG in exact arithmetic converges in one iteration. The cost analysis in Table 2 is consistent with dense tensor-mode products. The central risk is not the algebra but the empirical claim about iteration counts: the paper never defines the residual threshold that 'two or three sufficient' refers to, and it never reports iteration counts for the nine test problems. Figure 9 plots time, not iterations, and starts after initialization. This makes the headline claim unverifiable. A fixed-tolerance rerun would settle it. The reader's weakest assumption (storage of dense objects) is a real but explicitly acknowledged limitation (Section 6), and the reader's rationale already notes the post hoc iteration tuning; my concern partially overlaps with the reader's rationale but is not the reader's chosen weakest assumption. The verdict CONDITIONAL remains appropriate, and no adjustment is needed.","tokens_in":29308,"tokens_out":24930,"duration_ms":221668,"concrete_test":"Re-run all nine test cases from Problems 5.1–5.3 with the pseudoinverse preconditioner and a fixed stopping criterion, e.g., relative true residual ||H - L U_s||_F / ||H||_F ≤ 1e-10, and record the iteration count for each case. If the 500×1000 2D case or the 512×256×8 3D case requires more than 5 iterations to reach this tolerance, the Section 6 claim that 'two or three are sufficient' fails for the stated test set.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central efficiency claim — 'only a few iterations necessary ... two or three are sufficient' (Section 6) — is not backed by a quantitative stopping rule. Remark 4.2 says the iteration number is fixed experimentally 'so that the true residual decreases sufficiently,' but no target residual is specified. The experiments normalize ||R0|| to 1/N (Section 5.2), so for the largest 3D problem (512×256×8, N=1,048,576) the initial residual norm is ~1e-6; a vague 'sufficient' threshold could be orders of magnitude looser than what a typical solver user would require. Furthermore, Figure 9 reports wall-clock time (log scale) rather than iteration counts, and the curves begin after initialization, so the actual number of iterations for each of the nine test cases is never shown. The claim that two or three iterations are sufficient is therefore unfalsifiable from the presented data. Because the pseudoinverse preconditioner is mathematically the exact inverse on the range, a one-iteration convergence is expected in exact arithmetic; the practical question is how many iterations are needed in floating point at a specified tolerance, and the paper does not answer this.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper studies the preconditioned conjugate gradient (PCG) solution of finite-difference discretizations of Poisson's equation on rectangular or cuboid domains, where the discrete Laplacian is a sum of Kronecker products. The authors derive three preconditioners adapted to this structure: a Jacobi-like fixed-point preconditioner, a low-Kronecker-rank approximate inverse, and a Moore-Penrose pseudoinverse of the discretized Laplacian. The pseudoinverse preconditioner is built from analytic spectral decompositions of the one-dimensional Laplacians and is applied through basis transformations and Hadamard multiplication. Numerical experiments on six 2D and three 3D right-hand sides show rapid residual decrease, and the paper concludes that two or three PCG iterations are sufficient for this problem.","tokens_in":29538,"tokens_out":6593,"duration_ms":67288,"significance":"The mathematical derivation of the pseudoinverse preconditioner from the eigendecompositions of the 1D Laplacians is correct and standard, and the paper honestly identifies the method as a pseudoinverse-based direct solver with CG refinement. If the performance claims are substantiated, the method is a practical and simple Fourier-like preconditioner for Kronecker-structured Poisson problems with dense right-hand sides, and the connection to the Ferrodo2 application in ferroelectric simulation is valuable. The paper also provides explicit algorithms, analytical eigenvector formulas, and cost tables, which are useful for reproduction. However, the central efficiency claim is currently supported mainly by wall-clock plots with post hoc iteration counts and no comparison against standard baselines, so the significance of the numerical evidence is not yet at the level of the paper's conclusions.","major_comments":[{"comment":"The central claim in Section 6 that two or three PCG iterations are sufficient is not supported by the reported data. Remark 4.2 states that the iteration number is fixed experimentally so that the true residual decreases sufficiently, but no target residual or achieved residual is reported. Section 5.2 selects iteration counts post hoc, and Figure 9 shows only wall-clock time after initialization, with no iteration counts for any of the nine test cases. Since the pseudoinverse preconditioner is the exact inverse on the range of the operator, one iteration is expected in exact arithmetic; the practical question is how many iterations are needed in floating-point arithmetic at a specified relative tolerance. The paper does not answer this question, and the claim is therefore not reproducible from the presented data. Please provide a table of iteration counts and true residual norms relative to the normalized initial residual for each test problem.","section":"§5.2, Remark 4.2, Fig. 9"},{"comment":"The performance comparison lacks essential baselines. Figure 6 compares preconditioners only on a single 50×100 problem, and the larger 2D and 3D runs in Figure 9 report only pseudoinverse-preconditioned PCG times. No comparison is made with unpreconditioned sparse CG, a sparse direct solver, or an FFT-based direct solver, even though Remark 4.3 identifies the pseudoinverse preconditioner as essentially a Fourier solver. Without such baselines, the conclusion that the method is extremely efficient for practical problem sizes is not established. Please report timings, iteration counts, and achieved residuals against at least one standard sparse solver and one FFT/direct solver for the nine test problems.","section":"§5.2, Figs. 4–9"},{"comment":"The residual metric used in the experiments is not adequate for the application claim. Equation (5.1) normalizes the initial residual norm to 1/(nq) or 1/(nqt), so plots of absolute residual norms are not comparable across grids, and no relative tolerance is specified. In addition, the physically relevant output for the ferroelectric application is the electric field, i.e., the gradient of the potential, but the paper reports only residual norms and does not report error in the gradient. Please report relative residual norms and, where possible, the error in the computed electric field relative to a reference solution.","section":"§5.2, Eq. (5.1)"},{"comment":"The applicability of the method depends critically on the ability to store dense right-hand sides and dense intermediate matrices of size O(nq) or O(nqt). Section 6 acknowledges this, but the limitation should be stated more prominently, and the numerical experiments should include a memory comparison with sparse iterative solvers for the reported grid sizes. In regimes where dense storage is infeasible, the method loses its claimed advantage; the paper should make the scope of the efficiency claim explicit and quantify at which grid sizes the dense approach becomes impractical.","section":"§6"}],"minor_comments":[{"comment":"The phrase 'Newton condition e_R = -1/2' should be 'Neumann condition'.","section":"Problem 5.2"},{"comment":"The two formulas for the eigenvectors v_k are not explicitly assigned to L_{n,DN} and L_{n,ND}; please label them clearly.","section":"Lemma 3.4"},{"comment":"The indices ζ and ξ are used for entries of W without being defined; please replace them with i and j or define them in the notation section.","section":"§4.2"},{"comment":"There are small language issues: 'spaned' should be 'spanned', and 'Until stopping criterion match' should be 'Until stopping criterion matches'.","section":"§4.1 and Algorithm 1"},{"comment":"The paper states that no algorithmic description of the low-Kronecker-rank preconditioner is provided, yet this preconditioner is used in the experiments reported in Figures 6–8. Please either provide the algorithm or describe more precisely how the experiments were implemented, so that the comparison is reproducible.","section":"§4.3.2"},{"comment":"Please proofread for language issues such as 'preconditioners that fits', 'homogenous', and 'isotropous'.","section":"Abstract and Section 1"},{"comment":"The caption claims the convergence curve is 'almost problem independent', but the plot uses wall-clock time after initialization on a log scale and shows no iteration counts; please rephrase or provide the actual iteration counts and residual values.","section":"Fig. 9 caption"}],"recommendation":"major_revision","confidential_remarks":"The derivation of the pseudoinverse preconditioner is sound and the paper is honest about its dense-storage limitation. The main weakness is the numerical evaluation: post hoc iteration counts, no defined stopping tolerance, no FFT/direct-solver baseline, and no iteration counts for the larger problems. These issues are fixable within the scope of the manuscript, so I recommend major revision rather than rejection. I did not find citation or novelty concerns; the connection to the authors' PRB 2023 work is appropriate."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Dear [Colleague],\n\nThis is a case study, not a new algorithm, and the authors know it. The genuinely useful part is a Moore–Penrose pseudoinverse preconditioner for the Kronecker-structured finite-difference Laplacian on cuboids with periodic, Dirichlet, Neumann, and mixed boundary conditions, constructed from the analytic 1D spectra, plus cost tables. The derivation is correct and the exposition is clear. They also earn credit for Remark 4.3, which says plainly that the pseudoinverse preconditioner is essentially a Fourier solver—that is the right way to see it, and it undercuts any claim of novelty.\n\nThe central efficiency claim—\"two or three [iterations] are sufficient\"—is plausible in exact arithmetic, since the pseudoinverse is the inverse on the range. But the paper does not demonstrate it cleanly. Remark 4.2 admits the iteration count is fixed experimentally, with no target residual. Figure 9 reports wall-clock time starting after initialization, and iteration counts are never shown anywhere. So \"two or three\" is neither supported nor falsifiable from the data. The stress-test concern is correct.\n\nThe bigger soft spot is the missing baseline. For this separable Poisson problem on a rectangle, an FFT/DCT solver is the natural comparison, and per Remark 4.3 the method is a dense-transform Fourier solver. The cost analysis shows per-iteration work O(nq(n+q)) in 2D, which is worse than FFT asymptotically, but without measured constants or a comparison we cannot tell whether the mixed-boundary flexibility makes it competitive in the target application. The paper also requires dense storage of the n×q or n×q×t right-hand side and intermediate arrays, a limitation the authors state explicitly.\n\nNo code or data are released. For a methods case study, that is a real handicap: the runtime comparisons cannot be reproduced or checked.\n\nWho this is for: people building solvers for ferroelectric domain simulations, or those interested in Kronecker-structured preconditioning for a specific application. It deserves a serious referee, but the referee should require a real benchmark: iteration counts for every test case, a fixed stopping tolerance, and an FFT/DCT or direct-solver baseline. With those additions, the paper would become a useful, citable reference for its community.\n\nMy recommendation: send it to peer review, expecting major revisions on the experimental section.","headline":"A sound, honest case study whose central speed claim is plausible but not demonstrated by the reported experiments.","tokens_in":30086,"tokens_out":3045,"would_cite":false,"duration_ms":31127,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["15A06","15A30","15A69","65F08","65F10","65N06"],"pacs":[],"model":"deepseek-v4-flash","headline":"A pseudoinverse preconditioner built from the Kronecker structure of the finite-difference Laplacian lets PCG solve cuboid-grid Poisson problems in about two or three iterations.","keywords":["conjugate gradients","preconditioner","tensor","Kronecker product","Laplace operator","Poisson equation","ferroelectric domain structure","Moore-Penrose pseudoinverse"],"falsifier":"Take a dense, null-space-orthogonal right-hand side on a cuboid grid with periodic boundary conditions in all directions and run PCG with the pseudoinverse preconditioner; if the true residual has not dropped to near machine precision relative to the initial residual within three iterations (the paper's stated 'two or three are sufficient'), the central convergence claim is refuted. A concrete target: on the $512\\times 256\\times 8$ problem described in Section 5, count iterations to reach $\\|H-\\mathcal{L}U_s\\|_F/\\|H\\|_F \\le 10^{-10}$ and check whether the count is at most three.","tokens_in":29109,"feed_emoji":"⚡","tokens_out":8663,"duration_ms":72252,"temperature":0.7,"pith_summary":"This paper seeks to establish that preconditioned conjugate gradients (PCG) can solve the finite-difference Poisson equation on a rectangle or cuboid in essentially two or three iterations, provided the preconditioner is the Moore-Penrose pseudoinverse of the discrete Laplacian. Because the Laplacian separates into a sum of Kronecker products of small one-dimensional difference matrices, its full eigendecomposition is available analytically, and the pseudoinverse acts entrywise in that eigenbasis. The paper argues that this makes the method effectively a direct solver whose residual is polished by the outer CG iteration, at a per-iteration cost comparable to one multiplication by the Laplacian. It supports the claim with 2D and 3D experiments up to roughly a million unknowns, including singular periodic and Neumann cases where the pseudoinverse with thresholding handles the null space.","feed_headline":"Pseudoinverse preconditioner makes PCG converge in 2-3 iterations","feed_subtitle":"A Kronecker-structured Laplacian turns CG into an almost-direct solver on cuboid grids.","key_machinery":"The load-bearing object is the family of 1D finite-difference Laplacians $L_n(\\alpha,\\beta,\\gamma)$, whose four boundary-condition variants (periodic, Dirichlet, Neumann, mixed) have closed-form sine-and-cosine eigen-decompositions (Lemmas 3.1-3.4). The 2D/3D discrete Laplacian is the Kronecker sum of these 1D matrices, so its eigenvalues are sums of 1D eigenvalues and its eigenvectors are Kronecker products of 1D eigenvectors. The Moore-Penrose pseudoinverse preconditioner $M_{\\mathrm{p-inv}}$ exploits this by applying $V_n^T R V_q$, multiplying entrywise by the Hadamard pseudoinverse of the matrix of eigenvalue sums $(\\lambda_i+\\lambda_j)^{\\dagger}$, and transforming back with $V_n F V_q^T$ (Algorithms 4-7), which the paper identifies as a Fourier-solver step for the residual.","core_discovery":"The central discovery is that the Moore-Penrose pseudoinverse of the Kronecker-structured finite-difference Laplacian is an almost perfect preconditioner for PCG on this problem. Writing the 2D operator as $(I_q\\otimes L_n)+(L_q\\otimes I_n)$, the eigenvectors are Kronecker products of the 1D Laplacian eigenvectors and the eigenvalues are sums $\\lambda_i(L_n)+\\lambda_j(L_q)$; the pseudoinverse preconditioner transforms the residual into this eigenbasis, divides by the nonzero eigenvalue sums (treating sums below $10^{-13}$ as zero), and transforms back. In the authors' experiments the true residual is essentially at the desired level after the first iteration and fully converged after two or three, so the paper proposes viewing the scheme as a pseudoinverse-based direct solver with refinement by CG iteration. The same construction extends to 3D with a three-way tensor and triples of eigenvalue sums.","pith_inferences":["Because the pseudoinverse preconditioner is a trigonometric (Fourier-type) solver, on strongly anisotropic grids where one dimension is much smaller than the others the per-iteration cost $O(nqt(n+q+t))$ can exceed that of an FFT-based Poisson solver with $O(N\\log N)$ work; the paper's method is best suited to near-cubic grids where the dense-array storage is already affordable.","The same pseudoinverse-of-the-Kronecker-sum construction should transfer to other separable operators on cuboids, for example anisotropic diffusion with different coefficients per axis or Schrodinger operators with separable potentials, whenever the 1D factors have known eigendecompositions; the paper's experiments only cover the scalar homogeneous Laplacian.","The $10^{-13}$ threshold means the method is not a true direct solver for systems with very small nonzero eigenvalues; the number of CG iterations would grow if the right-hand side had significant energy in near-null directions, so one could probe robustness by constructing such right-hand sides explicitly."],"forward_implications":["For the simple Laplace operator on a cuboid, PCG with the pseudoinverse preconditioner needs only two or three iterations to converge; the authors run ten to be safe.","The method behaves as a direct solver: the pseudoinverse application produces an almost-exact solution and the CG outer loop removes the residual left by the eigenvalue-sum threshold.","The per-iteration cost is $O(nq(n+q))$ in 2D and $O(nqt(n+q+t))$ in 3D, comparable to the unpreconditioned PCG iteration itself, with initialization cost dominated by the eigendecompositions of the 1D Laplacians.","For singular Laplacians (periodic or Neumann in all directions), centering the data and thresholding the pseudoinverse keeps the iteration in the null-space-orthogonal subspace, so the solution is recovered up to the expected additive constant.","The approach has been implemented in C/C++ and used to produce published results on ferroelectric domain walls cited in the paper."],"supporting_citations":[{"why":"Supplies the closed-form sine/cosine eigen-decompositions of the four 1D Laplacian variants used to build the preconditioner.","marker":"[10]"},{"why":"Source for the 1D finite-difference matrices and the Kronecker-product identities that give the 2D/3D operators their tensor structure.","marker":"[6]"},{"why":"Provides the vectorization identities that recast the matrix/tensor equations as Kronecker-structured linear systems.","marker":"[27]"},{"why":"Defines the Moore-Penrose pseudoinverse and its properties, the basis of the central preconditioner.","marker":"[1]"},{"why":"Original CG method that the preconditioned variant accelerates.","marker":"[9]"},{"why":"Gives the Krylov-subspace context and convergence understanding used to justify the PCG framework.","marker":"[13]"},{"why":"Supplies the eigendecomposition machinery and operation counts behind the claimed computational cost.","marker":"[7]"},{"why":"Reports the physical ferroelectric domain-wall results computed by the implemented solver, motivating the method.","marker":"[15]"}],"fun_headline_variants":["Pseudoinverse PCG: near-direct solver for ferroelectric patterns","Kronecker preconditioner turns PCG into a 3-iteration solver","PCG with pseudoinverse: solve ferroelectric Poisson in 2 steps","Tensor PCG: pseudoinverse makes it an almost-direct solver"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The claim rests on being able to store and multiply dense $n\\times q$ (or $n\\times q\\times t$) arrays, because the right-hand side and all intermediate objects in the preconditioner are dense; if only sparse matrix storage fits in memory, this preconditioner's advantage over sparse iterative solvers disappears.","fun_headline_variants_meta":{"raw":{"variants":["Pseudoinverse PCG: near-direct solver for ferroelectric patterns","Kronecker preconditioner turns PCG into a 3-iteration solver","PCG with pseudoinverse: solve ferroelectric Poisson in 2 steps","Tensor PCG: pseudoinverse makes it an almost-direct solver"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000241,"raw_usage":{"total_tokens":1527,"prompt_tokens":956,"completion_tokens":571,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":572,"completion_tokens_details":{"reasoning_tokens":490}},"tokens_in":572,"tokens_out":571,"duration_ms":89558,"temperature":1.0,"reasoning_tokens":490,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T21:52:22.259268+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a dense, null-space-orthogonal right-hand side on a cuboid grid with periodic boundary conditions in all directions and run PCG with the pseudoinverse preconditioner; if the true residual has not dropped to near machine precision relative to the initial residual within three iterations (the paper's stated 'two or three are sufficient'), the central convergence claim is refuted. A concrete target: on the $512\\times 256\\times 8$ problem described in Section 5, count iterations to reach $\\|H-\\mathcal{L}U_s\\|_F/\\|H\\|_F \\le 10^{-10}$ and check whether the count is at most three.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Original CG method that the preconditioned variant accelerates."},{"cited_title":"Kol´ arov´ a: Spectral properties of discretized Laplace operator (in Czech), Diploma the- sis, Technical University of Liberec, 2024","cited_arxiv_id":null,"evidence_quote":"Supplies the closed-form sine/cosine eigen-decompositions of the four 1D Laplacian variants used to build the preconditioner."},{"cited_title":"Fiedler: Special matrices and their applications in numerical mathematics (2nd ed.), Dover Publications, Dover, 2008","cited_arxiv_id":null,"evidence_quote":"Source for the 1D finite-difference matrices and the Kronecker-product identities that give the 2D/3D operators their tensor structure."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the vectorization identities that recast the matrix/tensor equations as Kronecker-structured linear systems."},{"cited_title":"Ben-Israel, T","cited_arxiv_id":null,"evidence_quote":"Defines the Moore-Penrose pseudoinverse and its properties, the basis of the central preconditioner."},{"cited_title":"Liesen, Z","cited_arxiv_id":null,"evidence_quote":"Gives the Krylov-subspace context and convergence understanding used to justify the PCG framework."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the eigendecomposition machinery and operation counts behind the claimed computational cost."},{"cited_title":"Marton, M","cited_arxiv_id":null,"evidence_quote":"Reports the physical ferroelectric domain-wall results computed by the implemented solver, motivating the method."}],"review_version":1}