Riemannian Newton Methods with Statistical Inference

Xiaotian Chang; Qirui Hu; Yangdi Jiang; Nan Qiao

stat.MLOct 31, 2026

Riemannian Newton Methods with Statistical Inference

Authors: Xiaotian Chang, Qirui Hu, Yangdi Jiang, Nan Qiao

Organizations: Division of Mathematical Sciences, Nanyang Technological University. · School of Statistics and Data Science, Shanghai University of Finance and Economics. · School of Mathematics, Statistics and Mechanics, Beijing University of Technology, China.

Abstract

We develop online statistical inference for Polyak-Ruppert averages of inexact Newton iterates on Riemannian manifolds. A local asymptotic representation separates observation noise, solver randomization, and accumulated implementation bias, and yields one-pass covariance estimators and asymptotically valid Wald regions under local regularity and centering conditions. Independent fresh solves add a positive semidefinite term to the score-sandwich covariance; averaging solver replicas reduces this term in inverse proportion to their number. Slowly varying predictable gains instead cancel from the first order law. Product torus simulations show how Newton scaling shortens an ill-conditioned transient and how retaining solver variation restores coverage. Intrinsic means and principal subspaces provide geometric applications, with a streaming data illustration.

Figures & tables

		RSGD		Proposed, $m=4$		Exact Newton
$\kappa$	$n$	Cov. (%)	RMSE	Cov. (%)	RMSE	Cov. (%)	RMSE
1	256	95.8	18.99	95.2	21.00	95.8	18.99
1	1024	95.2	9.84	94.6	11.11	95.2	9.84
1	4096	96.2	4.63	94.4	5.38	96.2	4.63
30	256	75.4	29.60	96.4	21.60	96.6	18.85
30	1024	96.0	9.56	95.0	10.77	94.8	9.56

Table 1: Product torus comparison across spectral ratios and sample sizes. Coverage is for joint

95\%

Wald regions; geodesic RMSE is multiplied by

10^{3}

. Proposed uses

m=4

independent coordinate solves. Every entry uses 500 paired replications; solver budgets and checkpoint errors are in the appendix.

$\kappa$	Method	25%	50%	100%
1	RSGD	9.35	6.57	4.63
	Exact Newton	9.35	6.57	4.63
	Proposed, $m=4$	10.86	7.57	5.38
30	RSGD	9.20	6.67	4.67
	Exact Newton	8.80	6.51	4.66
	Proposed, $m=4$	11.01	7.99	5.53

Table 2: Geodesic RMSE along the

n=4096

averaging window, multiplied by

10^{3}

. The checkpoints contain 896, 1792 and 3584 retained observations after a burn-in of 512. Each row follows the same 500 paths with a fixed step constant; Proposed uses

m=4

		Coverage (%), $m=1$		RMSE $\times 10^{3}$
$\kappa$	$n$	Direction	Score only	$m=1$	$m=4$	$m=16$
1	256	95.8	76.2	27.01	21.00	19.57
1	1024	93.0	75.2	14.01	11.11	10.16
1	4096	95.6	80.0	6.61	5.38	4.85
30	256	95.4	78.4	26.94	21.60	19.21
30	1024	95.4	78.6	13.40	10.77	9.86

Table 3: Solver covariance correction and replication. The two coverage columns use the same proposed estimator with

m=1

: direction covariance versus the score only counterfactual that omits solver variation. Geodesic RMSEs for

m\in\{1,4,16\}

use 500 replications and are multiplied by

10^{3}

Appendix figures & tables5 assets

Supplementary material from the paper’s appendix.

Appendix

			Coverage (%)		RMSE $\times 10^{3}$
$n$	Method	$c$	Cov.	Score	25%	50%	100%
256	RSGD	3.2	95.8	95.8	36.35	26.29	18.99
256	Exact Newton	3.2	95.8	95.8	36.35	26.29	18.99
256	Proposed, $m=1$	0.8	95.8	76.2	46.82	36.69	27.01
256	Proposed, $m=4$	1.6	95.2	91.6	39.29	29.49	21.00
256	Proposed, $m=16$	1.6	95.8	94.4	36.85	26.79	19.57

Appendix

Table 4: Complete product torus results at

\kappa=1

. Cov. is implemented direction coverage for proposed and score-sandwich coverage for RSGD and exact Newton. Score is the score only counterfactual for the same estimate. All RMSEs are multiplied by

10^{3}

; fractions refer to the retained averaging window. The constant

c

was selected before evaluation. Each row has 500 replications.

			Coverage (%)		RMSE $\times 10^{3}$
$n$	Method	$c$	Cov.	Score	25%	50%	100%
256	RSGD	1.6	75.4	75.4	87.75	52.02	29.60
256	Exact Newton	6.4	96.6	96.6	36.18	26.87	18.85
256	Proposed, $m=1$	0.8	95.4	78.4	46.75	36.03	26.94
256	Proposed, $m=4$	6.4	96.4	93.6	42.95	31.40	21.60
256	Proposed, $m=16$	6.4	98.4	97.4	37.41	27.80	19.21

Appendix

Table 5: Complete product torus results at

\kappa=30

. Cov. is implemented direction coverage for proposed and score-sandwich coverage for RSGD and exact Newton. Score is the score only counterfactual for the same estimate. All RMSEs are multiplied by

10^{3}

; fractions refer to the retained averaging window. The constant

c

was selected before evaluation. Each row has 500 replications.

			Coverage (%)		RMSE $\times 10^{3}$
$n$	Method	$c$	Cov.	Score	25%	50%	100%
256	RSGD	0.2	0.0	0.0	450.73	412.01	347.57
256	Exact Newton	3.2	95.8	95.8	35.90	25.89	18.88
256	Proposed, $m=1$	0.8	94.4	79.0	49.39	37.28	27.19
256	Proposed, $m=4$	0.8	95.4	92.6	36.57	27.67	20.68
256	Proposed, $m=16$	3.2	96.6	95.4	37.68	26.91	19.29

Appendix

Table 6: Complete product torus results at

\kappa=300

. Cov. is implemented direction coverage for proposed and score-sandwich coverage for RSGD and exact Newton. Score is the score only counterfactual for the same estimate. All RMSEs are multiplied by

10^{3}

; fractions refer to the retained averaging window. The constant

c

was selected before evaluation. Each row has 500 replications.

$\kappa$	$n$	RSGD	Exact Newton	Proposed, $m=1$	Proposed, $m=4$	Proposed, $m=16$
1	256	0.045	0.046	0.068	0.102	0.205
1	1024	0.182	0.180	0.257	0.407	0.799
1	4096	0.711	0.719	1.023	1.581	3.191
30	256	0.045	0.046	0.066	0.105	0.207
30	1024	0.176	0.177	0.262	0.403	0.803
30	4096	0.692	0.708	1.020	1.590	3.181

Appendix

Table 7: Product torus median serial time per complete path, in milliseconds. Medians are over 30 randomized-order blocks of 100 paths per cell. Timing includes observation generation and uses aggregated coordinate sketches.

Feature group	Projector mass	95% interval
Size (radius/perimeter/area)	0.392	[0.254, 0.530]
Boundary irregularity	0.455	[0.390, 0.521]
Texture	0.023	[0.005, 0.041]
Smoothness and symmetry	0.092	[0.041, 0.143]
Fractal dimension	0.038	[0.004, 0.071]

Appendix

Table 8: WDBC leading-projector mass by morphology-based feature group. Intervals are pointwise 95% delta-method intervals under the conditional working model.

Covariance Estimation

Covariance estimation develops methods for learning dependence and uncertainty structure from noisy, incomplete, or high-dimensional data.

See all 66 papers

Bootstrap Inference for Riemannian Optimization Problems
Oct 31, 2026Xiaotian Chang, Qirui Hu, Yangdi Jiang +1Confidence Region Estimation Riemannian Optimization
Division of Mathematical Sciences, Nanyang Technological University. · School of Statistics and Data Science, Shanghai University of Finance and Economics. · School of Mathematics, Statistics and Mechanics, Beijing University of Technology, China.
Online Covariance Estimation for Averaged Stochastic Gradient Descent on Riemannian Manifolds
Oct 24, 2026Jingrong Guan, Liwei Jiang, Sen NaCovariance Estimation Riemannian Optimization
Department of Industrial Engineering and Management Sciences, Northwestern University · Edwardson School of Industrial Engineering, Purdue University · School of Industrial and Systems Engineering, Georgia Institute of Technology
Locally Private Inference for Riemannian Stochastic Optimization
Sep 18, 2026Xiaotian Chang, Yangdi Jiang, Qirui HuRiemannian Optimization Confidence Region Estimation
School of Physical and Mathematical Sciences, Nanyang Technological University · School of Statistics and Data Science, Shanghai University of Finance and Economics.

Oct 31, 2026J/K move · Enter open · S save

Bootstrap Inference for Riemannian Optimization Problems

Xiaotian Chang, Qirui Hu, Yangdi Jiang, Nan Qiao

Division of Mathematical Sciences, Nanyang Technological University. · School of Statistics and Data Science, Shanghai University of Finance and Economics. · School of Mathematics, Statistics and Mechanics, Beijing University of Technology, China.

Confidence Region Estimation Riemannian Optimization Stochastic Approximation

Open paper AskPDF arXiv

We develop an online multiplier bootstrap for manifold-valued parameters estimated by averaging the iterates of Riemannian stochastic approximation. The method runs weighted iterations alongside the original algorithm, with each path using the same observations while updating its own state and intrinsic average. The weighted averages provide uncertainty assessments without a plug-in asymptotic covariance estimator. Under local stability and regularity conditions, we establish bootstrap consistency and show that a fixed number of paths suffices for asymptotically valid Student intervals and intrinsic confidence regions. The theory accommodates samplewise nonsmooth scores and identifies a critical oracle-bias regime in which the estimator's limiting distribution shifts, but the centered bootstrap distribution does not. Applications include spherical means, positive-definite matrix means and medians, and principal-subspace estimation. Simulations show coverage close to nominal in regular settings, and an analysis of North Atlantic storm locations illustrates descriptive online inference for a spherical center.

Abstract

Covariance Estimation

Explore similar work

Bootstrap Inference for Riemannian Optimization Problems

Online Covariance Estimation for Averaged Stochastic Gradient Descent on Riemannian Manifolds

Locally Private Inference for Riemannian Stochastic Optimization