International Symposium on Recent Advances in Theories andMethodologies for Large Complex Data
December 7-9, 2023
Venue:
Conference Room 101, Tsukuba International Congress Center
2-20-3 Takezono, Tsukuba, Ibaraki 305-0032, Japan (Hybrid Symposium with Zoom)
Organizers:
Makoto Aoshima (University of Tsukuba)
Kazuyoshi Yata

(University of Tsukuba)

Aki Ishii(Tokyo University of Science)

Kento Egashira(Tokyo University of Science)

Supported by
Grant-in-Aid for Scientiﬁc Research (A) 20H00576 (Project Period: 2020-2024)“Innovative developments of theories and methodologies for large complex data”(Principal Investigator: Makoto Aoshima)
Grant-in-Aid for Challenging Research (Exploratory) 22K19769 (Project Period: 2022-2024)“Developments of statistical compression technology for massive data having tensor structures”(Principal Investigator: Makoto Aoshima)

Program (UTC+9)
December 7 (Thursday)
13:50∼14:00

Opening

14:00∼14:40

Kento Egashira∗,a , Kazuyoshi Yatab and Makoto Aoshimabab(Department of Information Sciences, Tokyo University of Science)(Institute of Mathematics, University of Tsukuba)

Asymptotic properties of kernel k-means for high dimensional data
14:50∼15:30

Yohji Akama(Mathematical Institute, Tohoku University)

Broken-stick components retention rule for equi-correlated normal population
15:40∼16:20

Takayuki Morimoto(School of Science, Kwansei Gakuin University)

Forecasting high-dimensional covariance matrices using high-dimensional
principal component analysis

16:30∼17:10

Shao-Hsuan Wang
(Graduate Institute of Statistics, National Central University)

A geometric algorithm for contrastive principal component analysis in high dimension
17:20∼18:00

Taiji Suzuki∗,a , Denny Wub , Atsushi Nitandac and Kazusato Okoa

(Zoom)abc(Department of Mathematical Informatics, The University of Tokyo / RIKEN AIP)(Center for Data Science, New York University)(Department of Artiﬁcial Intelligence, Kyushu Institute of Technology / RIKEN AIP)

Feature learning via mean ﬁeld neural networks and anisotropic features

December 8 (Friday)
9:00∼9:40

Akifumi Okuno(The Institute of Statistical Mathematics / RIKEN AIP)

Statistical estimation with integral-based loss functions
9:50∼10:30

Masaaki Imaizumi(Komaba Institute for Science, The University of Tokyo / RIKEN AIP)

Non-sparse high-dimensional statistics and its applications
10:40∼11:20

Tsutomu T. Takeuchi(Division of Particle and Astrophysical Science, Nagoya University)

Statistical challenges to dimensionality in astronomical big data
11:30∼12:10

Yuan-Tsung Chang∗,a , Nobuo Shinozakib and William, E. Strawdermancabc(Department of Social Information, Mejiro University)(Faculty of Science and Technology, Keio University)(Department of Statistics, Rutgers University)

Predictive density estimation for two ordered normal means under α-divergence loss
12:10∼13:40

Lunch

13:40∼18:00

Special Invited and Keynote Sessions

19:00∼21:00

DinnerDecember 9 (Saturday)
9:00∼9:40

Shogo Nakakita(Komaba Institute for Science, The University of Tokyo)

On approximate sampling from non-log-concave non-smooth distributions
via a Langevin-type Monte Carlo algorithm

9:50∼10:30

Kou Fujimori∗,a and Koji Tsukudabab(Department of Economics, Shinshu University)(Faculty of Mathematics, Kyushu University)

Two step estimations via the Dantzig selector for ergodic time series models
10:40∼11:20

Junichi Hirukawa∗,a and Kou Fujimoribab(Faculty of Science, Niigata University)(Department of Economics, Shinshu University)

Innovation algorithm of fractionally integrated (I(d)) process and applications
on the estimation of parameters
11:30∼12:10

Kengo Kamatani(The Institute of Statistical Mathematics)

Scaling limits of Markov chains/processes in Monte Carlo methods
12:20∼13:00

Takahiro Nishiyama∗,a and Masashi Hyodobab(Department of Business Administration, Senshu University)(Department of Economics, Kanagawa University)

On a general linear hypothesis testing problem for latent factor models
in high dimensions
13:00∼13:10

Closing(∗ Speaker)

Special Invited and Keynote Sessions
December 8 (Friday)
Special Invited Session
13:40∼14:30

On the eﬃciency-loss free ordering-robustness of product-PCA

Speaker: Hung Hung(Institute of Health Data Analytics and Statistics, National Taiwan University)
Discussion Leader: Yuan-Tsung Chang

14:40∼15:30(Department of Social Information, Mejiro University)

Learning ordinality in high-dimensional data

Speaker: Jeongyoun Ahn(Department of Industrial and Systems Engineering, KAIST)
Discussion Leader: Kazuyoshi Yata(Institute of Mathematics, University of Tsukuba)

Keynote Session
15:50∼16:50

Normal-reference test for high-dimensional covariance matrices

Speaker: Jin-Ting Zhang(Department of Statistics and Data Science, National University of Singapore)
Discussion Leader: Aki Ishii

17:00∼18:00(Department of Information Sciences, Tokyo University of Science)

Testing high-dimensional general linear hypotheses through spectral shrinkage

Speaker: Debashis Paul
(Zoom)(Department of Statistics, University of California, Davis / Indian Statistical Institute, Kolkata)

Discussion Leader: Yuta Koike(Graduate School of Mathematical Sciences, The University of Tokyo)

Asymptotic properties of kernel k-meansfor high dimensional data
Kento Egashiraa , Kazuyoshi Yatab , Makoto Aoshimaba

Department of Information Sciences, Tokyo University of SciencebInstitute of Mathematics, University of Tsukuba

Cluster analysis can be divided into two types: hierarchical and partitional. Hierarchical clustering groups data into dendrograms based on their cluster similarities
determined by a preset linkage function. A dendrogram enables the observation of
the process of merging or dividing clusters. For discussions on hierarchical cluster
analyses, see the works of Everitt et al. [4] and Hastie et al. [6], among others. Partitional clustering, as its name suggests, divides data into a pre-determined numberof clusters. K-means can be given on behalf of partitional clustering. Notably, kmeans has been approved as a useful tool for analyzing microarray gene expression
data. A characteristic of such data is that the number of variables was considerably larger than the sample size, giving high-dimensional, low-sample-size (HDLSS)
scenarios. Substantial work on HDLSS asymptotic clustering has been performedin recent years. For example, Liu et al. [8] proposed a two-way split statisticalsigniﬁcance-of-clustering (SigClust) method for HDLSS data. Ahn et al. [1] proposed hierarchical divisive clustering for high-dimensional asymptotics. Huang et al.
[5] modiﬁed SigClust using a soft thresholding approach. Kimes et al. [7] proposeda method for sequentially testing the statistical signiﬁcance of hierarchical clustering by controlling the family-wise error rate in HDLSS settings. Yata and Aoshima[10] presented the consistency properties of sample principal component scores andapplied them to clustering in high-dimensional settings. Nakayama et al. [9] investigated HDLSS clustering using kernel principal component analysis. Borysov et al.
[2] studied the behaviors of hierarchical clustering under several asymptotic settingsfrom a moderate dimension for HDLSS; however, the theoretical assumptions wereconsidered to be strict for HDLSS data owing to several simultaneous asymptotic
settings. Egashira et al. [3] explores practical assumptions to indicate the behavior of hierarchical clustering and obtained theoretical results in multiclass settings.
Given this background, asymptotic properties of k-means in the HDLSS settingsseems to have not been studied suﬃciently.
In this talk, we investigated k-means when both the dimension and sample size
approach inﬁnity at ﬁrst. Then, we explored kernel k-means in the HDLSS context
theoretically. Especially, we mensioned kernel k-means with gaussian kernel function
and compared performance of it to conventional k-means in the multiclass HDLSS
context.

References[1] Ahn, J., Lee, M.H., Yoon, Y.J. (2012). Clustering high dimension, low sample
size data using the maximal data piling distance. Statistica Sinica, 22, 443–464.
[2] Borysov, P., Hannig, J., Marron, J.S. (2014). Asymptotics of hierarchical clustering for growing dimension. Journal of Multivariate Analysis, 124, 465–479.
[3] Egashira, K., Yata, K., Aoshima, M. (2023). Asymptotic properties of hierarchical clustering in high-dimensional settings. Journal of Multivariate Analysis,
199, 105251.
[4] Everitt, B.S., Landau, S., Leese, M. (2001). Cluster Analysis. Arnold, NewYork.
[5] Huang, H., Liu, Y., Yuan, M., Marron, J.S. (2015). Statistical Signiﬁcance ofClustering using Soft Thresholding. Journal of Computational and Graphical
Statistics, 24, 975–993.
[6] Hastie, T., Tibshirani, R., Friedman, J. (2009). The Elements of Statistical
Learning: Data Mining, Inference, and Prediction, Second Edition, Springer,
New York.
[7] Kimes, P.K., Liu, Y., Neil, H.D., Marron, J.S. (2017) Statistical signiﬁcance for
hierarchical clustering. Biometrics, 73, 811–821.
[8] Liu, Y., Hayes, D.N., Nobel, A., Marron, J.S.(2008). Statistical signiﬁcance ofclustering for high-dimension, low-sample size data. Journal of the American
Statistical Association, 103, 1281–1293.
[9] Nakayama, Y., Yata, K., Aoshima, M. (2021). Clustering by principal component analysis with Gaussian kernel in high-dimension, low-sample-size settings.
Journal of Multivariate Analysis, 185, 104779.
[10] Yata, K., Aoshima, M. (2020). Geometric consistency of principal componentscores for high-dimensional mixture models and its application. Scandinavian
Journal of Statistics, 47, 899–921.

Broken-stick components retention rule for equi-correlated
normal population (報告書)赤間陽二

2023 年 12 月 9 日
株価の低頻度サンプリング時系列や, 大規模心理学的調査などでは, 変数の間の相関が顕著になる. そのような状況で因子モデルの因子の個数を標本相関行列 C から推定することを考える. 代表的な推定方法としては,
閾値法 (thresholding) の他に, 「ランダム」データから因子などをシミュレーションして, C の対応物と比較する方法がある. 閾値法としては, Guttman-Kaiser rule や adjusted correlation thresholding [6] があり, 後者のシミュレーションによる因子数推定法としては, 伝統的な broken-stick rule や, 頻繁に使用される parallel

analysis [3] がある. 閾値法もシミュレーションによる因子推定法も, 実データを用いて振る舞いが研究されてきた. 閾値法の理論的背景は [8], [6], [2] などが考察しているが, broken-stick rule や parallel analysis の理論
的背景は十分理解されてない [5].
そこで broken-stick rule の理論研究のために次のような設定を考えた. どの異なる変数の組みも共通の定数 ρ ∈ [0, 1) を相関係数として持つ p 次元正規母集団 (等相関正規母集団と呼ぶことにする) に着目
し, サイズ n のサンプルのサンプル相関行列 C に対する broken-stick rule の結果の, 比例的漸近枠組み

n, p → ∞, p/n → c > 0 での極限を研究することにした. このような C の絶対最大非対角要素は, さまざまな漸近枠組みで ρ に関する相転移を現象を持つ [7]. 我々は, ρ の 0 または正で, Guttman-Kaiser rule の極限
挙動に相転移 [2] を見たように, broken-stick rule の極限挙動に相転移を見ることを目標とした.
本研究では, 比例的漸近枠組みにおいて, 等相関正規母集団の標本相関行列 C の第二固有値が (1 −ρ)(1 +√c)2 に高確率で収束することを証明した. その証明には, 比例的漸近枠組みでの非有界スペクトルを持つ標本分散行列 S の固有値分布 [4], S と C の漸近的関係 [1], および, Weyl の不等式を用いた. この第二固有値に関する命題と, C の最大固有値の発散 [1] より, broken-stick rule が C から抽出する因子の個数が,性質 ρ > 0 の指示関数に高確率で収束することが証明できた. この定理により, broken-stick rule は, 母相関係数が正である equi-correlation block の個数を数えていることが説明できた. 対応しそうな現象が, binary

multiple sequence alignment のデータ [9] や, Fama-French 100 portfolios のデータ [6] の平均相関係数の低頻度サンプリング時系列に確認できた.
しかし, シンポジウムで高次元ファイナンスの専門家が次のような指摘をした: 高頻度サンプリング時系列
分析では, 正規分布より裾が重い分布, 例えば t-分布などを用い, 因子モデルの因子の個数を推定するのに [10]の Shrinkage Principal Orthogonal complEment Thresholding (S-POET) という閾値法が有用である.
さらに, シンポジウムで矢田・青嶋のグループは, 一般的なスパイク固有値モデルの標本分散行列の固有値に関する分布自由な結果や, 比例的漸近枠組み以外の漸近枠組みに関する結果を示唆し, S-POET は noise

reduction PCA [11, 12] に関係していると指摘した.
以上から Yata & Aoshima や Fan の結果が, 本研究の理論展開に重要であることを認識した. またシンポ1

ジウムでは本研究に有用なさまざまな知見を得た.

参考文献[1] Y. Akama. Correlation matrix of equi-correlated normal population: ﬂuctuation of the largest eigenvalue, scaling of the bulk eigenvalues, and stock market. Int. J. Theor. Appl. Finance, 26:2350006,
2023.
[2] Y. Akama and A. Husnaqilati. A dichotomous behavior of Guttman-Kaiser criterion from equicorrelated normal population. J. Indones. Math. Soc., 28(3):272–303, 2022.
[3] A. Buja and N. Eyuboglu. Remarks on parallel analysis. Multiv. Behav. Res., 27:509–540, 1992.
[4] T. Cai, X. Han, and G. Pan. Limiting laws for divergent spiked eigenvalues and largest nonspiked
eigenvalue of sample covariance matrices. Ann. Stat., 48(3):1255–1280, 2020.
[5] E. Dobriban and A. B. Owen. Deterministic Parallel Analysis: An Improved Method for Selecting
Factors and Principal Components. J. R. Stat. Soc. Ser. B Stat. Methodol., 81(1):163–183, 11 2018.
[6] J. Fan, J. Guo, and S. Zheng. Estimating number of factors by adjusted eigenvalues thresholding.
J. Am. Stat. Assoc., 117(538):852–861, 2022.
[7] J. Fan and T. Jiang. Largest entries of sample correlation matrices from equi-correlated normal
populations. Ann. Probab., 47(5):3321–3374, 2019.

√[8] M. Gavish and D. L. Donoho. The optimal hard threshold for singular values is 4/ 3. IEEE Trans.
Inf. Theory, 60(8):5040–5053, 2014.
[9] A. A. Quadeer, R. H. Louie, K. Shekhar, A. K. Chakraborty, I. Hsing, and M. R. McKay. Statistical linkage analysis of substitutions in patient-derived sequences of genotype 1a hepatitis C virus
nonstructural protein 3 exposes targets for immunogen design. J. Virol., 88(13):7628–7644, 2014.
[10] Weichen Wang and Jianqing Fan. Asymptotics of empirical eigenstructure for high dimensional
spiked covariance. The Annals of Statistics, 45(3):1342 – 1374, 2017.
[11] Kazuyoshi Yata and Makoto Aoshima. Eﬀective pca for high-dimension, low-sample-size data withnoise reduction via geometric representations. Journal of Multivariate Analysis, 105(1):193–215,
2012.
[12] Kazuyoshi Yata and Makoto Aoshima.

Pca consistency for the power spiked model in high-dimensional settings. Journal of Multivariate Analysis, 122:334–354, 2013.

2

Forecasting High-Dimensional CovarianceMatrices Using High-Dimensional Principal
Component Analysis
Hideto Shigemoto ∗, Takayuki Morimoto †November 17, 2023

Abstract
We modify the recently proposed forecasting model of highdimensional covariance matrices (HDCM) of asset returns using
high-dimensional principal component analysis (PCA). It is wellknown that when the sample size is smaller than the dimension,
eigenvalues estimated by classical PCA have a bias. In particular, avery small number of eigenvalues are extremely large and they arecalled spiked eigenvalues. High-dimensional PCA gives eigenvalues
which correct the biases of the spiked eigenvalues. This situation
also happens in the ﬁnancial ﬁeld, especially in situations where
high-frequency and high-dimensional data are handled. The researchaims to estimate the HDCM of asset returns using high-dimensional
PCA for the realized covariance matrix using the Nikkei 225 data,
it estimates 5- and 10-minute intraday asset-returns intervals. We
construct time-series models for eigenvalues which are estimated byeach PCA, and forecast HDCM. Our simulation analysis shows thatthe high-dimensional PCA has better estimation performance thanclassical PCA for the estimating integrated covariance matrix. In
our empirical analysis, we show that we will be able to improve theforecasting performance using the high-dimensional PCA and makea portfolio with smaller variance.

∗ Group

Risk Management Department, Nomura Holdings, Inc., 2-2-2 Otemachi,
Chiyoda-ku, Tokyo 100-8130, Japan.
† Corresponding author: Department of Mathematical Sciences, Kwansei Gakuin University, 1, Gakuen Uegahara, Sanda, Hyogo 669-1330, Japan.
morimot@kwansei.ac.jp

1Keywords: covariance forecasting; high-dimensional covariance; principal
component analysis; high-frequency data; time series

Summary
In this study, we constructed the HDCM forecasting models using highdimensional PCA. In particular, the previous studies show that to estimate
the latent factors, POET is used. However, it is known that when the dimension is greater than the sample size, the eigenvalues estimated by classical
PCA have biases. Therefore, in order to estimate the eigenvalues more accurately, we adopted SPOET which corrects biases of empirical eigenvalues. In
addition, we combined eigenvalues and time-series models to forecast eigenvalues and covariance matrix.
In the simulation study, we generated the asset returns based on the estimated HDCM as the integrated covariance matrix and it shows that SPOET
is also eﬀective for the price process. Especially, the empirical eigenvalues of
SPOET were closer to the true values than those of POET.
In the empirical analysis, we constructed some forecasting models of
HDCM using a number of individual stocks traded on Nikkei 225. Almost
all our proposed models which use SPOET show better performance thanthe other models which use POET. In addition, in terms of economic
performance, our models can generate a smaller variance than benchmarks
in most cases. This study applied SPOET discussed under the i.i.d. setting
to the continuous Itô semi-martingale setting for simulation study andempirical analysis. Thus, theoretical results are needed in the future.

This study is partly supported by the Institute of Statistical Mathematics(ISM) cooperative research program (2022-ISMCRP-2024), JSPS KAKENHI
Grant Number 21K01433, and Grant-in-Aid for JSPS Fellows Grant Number
22J10285. The views expressed in this study are those of the authors and do
not necessarily reﬂect the oﬃcial views of Nomura Holdings, Inc. or Kwansei
Gakuin University.

2$*HRPHWULF$OJRULWKPIRU&RQWUDVWLYH3ULQFLSDO&RPSRQHQW$QDO\VLV
LQ+LJK'LPHQVLRQ

6KDR+VXDQ:DQJ*UDGXDWH,QVWLWXWHRI6WDWLVWLFV1DWLRQDO&HQWUDO8QLYHUVLW\$EVWUDFW 3ULQFLSDO FRPSRQHQW DQDO\VLV 3&$  KDV EHHQ ZLGHO\
XVHGLQH[SORUDWRU\GDWDDQDO\VLV&RQWUDVWLYH3&$ $ELGHWDO DJHQHUDOL]HGPHWKRGRI3&$LVDQHZWRROXVHGWRFDSWXUHIHDWXUHV RI D WDUJHW GDWDVHW UHODWLYH WR D EDFNJURXQG GDWDVHWZKLOHSUHVHUYLQJWKHPD[LPXPDPRXQWRILQIRUPDWLRQFRQWDLQHGLQ
WKH GDWD :LWK KLJK GLPHQVLRQDO GDWD FRQWUDVWLYH 3&$ EHFRPHVLPSUDFWLFDOGXHWRLWVKLJKFRPSXWDWLRQDOUHTXLUHPHQWRIIRUPLQJWKH FRQWUDVWLYH FRYDULDQFH PDWUL[ DQG DVVRFLDWHG HLJHQYDOXHGHFRPSRVLWLRQIRUH[WUDFWLQJOHDGLQJFRPSRQHQWV,QWKLVZRUNZHSURSRVHDJHRPHWULFFXUYLOLQHDUVHDUFKPHWKRGWRVROYHWKLVSUREOHPDQGSURYLGHDFRQYHUJHQFHDQDO\VLV2XUDSSURDFKRIIHUVVLJQLILFDQWFRPSXWDWLRQDOHIILFLHQFLHV6SHFLILFDOO\LWUHGXFHVWKHWLPHFRPSOH[LW\IURP2 QҎP SA WRDPRUHPDQDJHDEOH2 Q
ҎP SU ZKHUHQPDUHWKHVDPSOHVL]HVRIWKHWDUJHWGDWDDQG
EDFNJURXQGGDWDUHVSHFWLYHO\SLVWKHGDWDGLPHQVLRQDQGULV
WKH QXPEHU RI OHDGLQJ FRPSRQHQWV $GGLWLRQDOO\ ZH VWUHDPOLQH
WKH VSDFH FRPSOH[LW\ IURP 2 SA  QHFHVVDU\ IRU VWRULQJ WKHFRQWUDVWLYHFRYDULDQFHPDWUL[WRDPRUHHFRQRPLFDO2 QҎP S VXIILFLHQW IRU VWRULQJ WKH GDWD DORQH 1XPHULFDO H[DPSOHV DUH
SUHVHQWHGWRVKRZWKHPHULWVRIWKHSURSRVHGDOJRULWKP

Feature learning via mean-ﬁeld neural networks andanisotropic features
Taiji Suzuki1,2
Joint work with Denny Wu3 , Atsushi Nitanda2,4 , Kazusato Oko1,2
1

Graduate School of Information Science and Technology, the University of Tokyo,
2
RIKEN Center for Advanced Intelligence Project,
3
Center for Data Science, New York University,
4
Faculty of Computer Science and Systems Engineering, Kyushu Institute of Technology.
In this presentation, I have presented the feature learning perspective of mean-ﬁeld neuralnetworks and compared it with kernel methods. Especially, I presented the generalization ability ofmean ﬁeld neural network trained by mean ﬁeld Langevin dynamics for learning k-sparse parity
functions. I summarize the technical details given in the presentation as follows.

1

Introduction

In this work, we bridge the aforementioned gap by presenting a simple and general framework toestablish sample complexity of MFLD in learning binary classiﬁcation problems. We then apply
this framework to the sparse k-parity problem, and obtain improved rate of convergence for the
fully time- and space-discretized algorithm. More speciﬁcally, our contributions can be summarized
as follows.
• We present a general framework to analyze MFLD in the learning of binary classiﬁcation tasks.
Our framework has two main ingredients: (i) an annealing procedure that applies to common
classiﬁcation losses that removes the exponential dependence on regularization parameters in
the logarithmic Sobolev inequality, and (ii) a novel local Rademacher complexity analysis for the
distribution of parameters optimized by MFLD. As a result, we can obtain learning guarantees
for the mean-ﬁeld neural network in discrete-time and ﬁnite-width settings.
• We apply our general framework to the k-sparse parity problem, and derived learning guarantees
with improved rate of convergence and dimension dependence. Specially, in the n  d2 regime
we obtain exponentially converging classiﬁcation error, whereas in the n  d regime we achieve
linear dimension dependence. Note that this improves upon the NTK analysis (which gives asample complexity of n = Ω(dk )) in that it “decouples” the degree k from the exponent in thedimension dependence. Our theoretical results are supported by empirical ﬁndings.

2

Problem Setting

We consider a classiﬁcation problem given by the following model:
Y = 1A (Z) − 1Ac (Z) ∈ {±1}where Z = (Z1 , . . . , Zd ) is the input random variable on Rd and 1A is the indicator function
corresponding to a measurable set A ∈ B(Rd ), i.e., 1A (Z) = 1 if Z ∈ A and 1A (Z) = 0 if Z ∈ A.
Let PZ be the distribution of Z. We are given input-output pairs Dn = (zi , yi )ni=1 independently
identically distributed from this model as training data. Then, we construct a binary classiﬁer
that predicts the label for the test input data as accurate as possible. To achieve this, we learn a
two-layer neural network model in the mean-ﬁeld regime via the mean-ﬁeld Langevin dynamics.
One important problem setting for our analysis is the k-sparse parity problem deﬁned as follows.

1

√
Example 1 (k-sparse parity √problem). PZ is the uniform distribution on the grid {±1/ d}d and
A = {ζ = (ζ1 , . . . , ζd ) ∈ {±1/ d}d | ζ1 · · · ζk > 0}1 .
Mean-ﬁeld two-layer network. Given input z, let hx (z) be one neuron in a two-layer neuralnetwork with parameter x = (x1 , x2 , x3 ) ∈ Rd+1+1 deﬁned ashx (z) = R̄[tanh(z  x1 + x2 ) + 2 tanh(x3 )]/3,
where R̄ ∈ R is a hyper-parameter determining the scale of the network. We place an extra tanh
activation for the bias term x3 ∈ R because the boundedness of hx is required in the convergence¯¯¯analysis. Let P be the set of probability measures on (Rd , B(Rd )) where d¯ = d + 2 and B(Rd )d¯is the Borel σ-algebra on R and Pp be the subset of P such that its p-th moment is bounded:Eμ [Xp ] < ∞ (μ ∈ P). The mean-ﬁeld neural network is deﬁned as an integral over neurons hx ,fμ (·) = hx (·)μ(dx),for μ ∈ P. To evaluate the performance of fμ , we deﬁne the empirical risk and the population riskasnL(μ) := n1 i=1 (yi fμ (zi )), L̄(μ) := E[(Y fμ (Z))],
respectively, where  : R → R≥0 is a convex loss function. In particular, we consider the logistic loss
(f, y) = log(1 + exp(−yf )) for y ∈ {±1} and f ∈ R. To avoid overﬁtting, we consider a regularizedempirical risk F (μ) := L(μ) + λEX∼μ [λ1 X2 ], where λ, λ1 ≥ 0 are regularization parameters.
One advantage of this mean-ﬁeld deﬁnition is that fμ is a linear with respect to μ, and hence thefunctional L(μ) becomes a convex functional. Let μ̂ be the minimizer of F (μ). We put the following
assumption.
Assumption 1. There exists c0 > 0 and R > 0 such that the following conditions are satisﬁed:• For some R̄, there exists μ∗ ∈ P such that KL(ν, μ∗ ) ≤ R and L(μ∗ ) ≤ (0) − c0 .
• For any λ < c0 /R, the regularized expected risk minimizer μ[λ] := argmin L̄(μ) + λKL(ν, μ)satisﬁes Y fμ[λ] (X) ≥ c0 almost surely.
Type I: Perfect Classiﬁcation with Exponentially Decaying Error Under the margin
assumption of fμ∗ (Assumption 1), we have that fμ̂ also yields a Bayes optimal classiﬁer. Moreprecisely, we have the following theorem.
Theorem 1. Suppose Assumption 1 holds. Let M0 = (∗ + 2(R̄ + 1))/λ. Moreover, suppose that
λ < c0 /R and
 


4R̄2
λ
2
2
2
∗
Q := c0 −λ 4R̄ +
+ 8R̄ (4 + log log2 (8n M0 R̄)) + nλ > 0,nλ2
32R̄2 n2nλthen fμ̂ yields perfect classiﬁcation, i.e., P (Y fμ̂ (Z) > 0) = 1, with probability 1 − exp(− 32Q).
R̄4Type II: Polynomial Order Classiﬁcation Error We can also obtain a result that has amilder dependency on λ and hence a better sample complexity.
Theorem 2. Suppose Assumption 1 holds. Let λ < c0 /R and M0 = (∗ + 2(R̄ + 1))/λ. Then, withprobability 1 − exp(−t), the classiﬁcation error of fμ̂ is bounded as
 2



8R̄ 
1
λ
∗
P (Y fμ̂ (Z) ≤ 0) ≤ 2ψ(c0 ).
4 + t + log log2 (8n2 M0 R̄) +
4R̄ +
+

nλ
n
32R̄2 n
We notice that the right hand side scales with O(1/(nλ)), which is better than O(1/(nλ2 )) in
Theorem 1; this implies that a sample size linear in the dimensionality is suﬃcient to achieve smallclassiﬁcation error. The reason for such improvement in the λ-dependence is that the stronger
L∞ -norm convergence is not used in the proof; instead, only the convergence of the loss is utilized.
On the other hand, this analysis does not guarantee a perfect classiﬁcation.
1 We present the axis-aligned setting for conciseness, but the same result holds under orthogonal transforms.

2

Statistical estimation with integral-based loss functions
Akifumi Okuno1,2
1 Inst. Stat. Math., 2 RIKEN AIPThis study considers statistical estimation problems using integral-based loss functions,which encompass the following two speciﬁc applications.
(1) Robust estimation using density-power divergence discussed in Okuno (2023a). Therein,
they minimize the following robust loss function equipped with a parametric density p θ arbitrarily speciﬁed by users (e.g., Gaussian mixture model, Gompertz model):

2

4

6

8

0.20
0.15
0.10
0.05
0.00

0

0.0

0
0

10

Probability density6080

MLE (beta=0)
beta=0.5 / m=10
beta=0.5 / m=25
beta=0.5 / m=10020

Frequency

0.25

p θ (z)1+β dz.

40

0.1

Probability density

0.4
0.3
0.2

60
40
20

Frequency80

MLE (beta=0)
beta=1 / m=10
beta=1 / m=25
beta=1 / m=100



100 120n11
1−p θ (x i )β +β n i =1
1+β

−5

0

5

10

x

x

(a) Gompertz density(b) Gaussian mixture

Figure 1: Computational intractable models (Gompertz and Gaussian mixture) are estimated by
the proposed approach.
(2) Higher-order variation regularization discussed in Okuno (2023b). Therein, they minimize
the following penalized loss function equipped with a function f θ arbitrarily speciﬁed byusers (e.g., neural network, generalized additive model):

−0.4

−0.2

0.0

0.2

0.4

x

(a) No regularization

1.0
y

0.0

0.2

0.4

0.6

0.8

1.0
0.8
0.6
y

0.4
0.0

0.2

0.4
0.0

0.2

y

0.6

0.8

1.0

 k
q
n
K

1
∂

2
{y i − f θ (x i )} +
η k  k f θ (z) dz.
n i =1
∂z
k=0

−0.4

−0.2

0.0

0.2

0.4

x

(b) Ridge regularization

−0.4

−0.2

0.0

0.2

0.4

x

(c) Proposal

Figure 2: Single-layer perceptron (with L = 200 hidden units) trained with the proposed approach.

Due to the difﬁculty in evaluating the integral terms, much of the existing research hasconcentrated on (i) straightforward models (e.g., normal density estimation, spline estimation)
where the explicit form of these integral terms can be obtained, or (ii) numerical integration,
which is computationally intensive. Notably, statisticians have long focused on optimization
with full-batch methods (e.g., Newton-Raphson method), and the obsession with full-batch methods has made computations challenging.
However, if we rewind the long history of research all the way back to the beginning, suchintegral-based loss functions have been known to be simply minimized by stochastic gradientdescent (Robbins and Monro, 1951). To bridge the substantial, unseen gap between the practicalapplications in statistics and the profound advancements in stochastic optimization theory, thistalk intentionally sheds light on the potential utility of the stochastic optimization techniques.

References
Okuno, A. (2023a). A stochastic optimization approach to minimize robust density power-based
divergences for general parametric density models. arXiv preprint arXiv:2307.05251.
Okuno, A. (2023b). A stochastic optimization approach to train non-linear neural networks witha higher-order variation regularization. arXiv preprint arXiv:2308.02293.
Robbins, H. and Monro, S. (1951). A Stochastic Approximation Method. The Annals of Mathematical Statistics, 22(3):400 – 407.

NON-SPARSE HIGH-DIMENSIONAL STATISTICS AND ITS APPLICATIONSMASAAKI IMAIZUMI1,2
1 The University of Tokyo, 2 RIKEN Advanced Intelligence Project

Abstract. In this talk, we present several results in non-sparse high-dimensional statistics. Specifically, the generalization and Bayesian estimation of high-dimensional linear regression models,statistical inference for high-dimensional generalized linear models, and the regret analysis on contextual bandit problems applying high-dimensional linear models. The analysis in these studies uses
the theory of benign overﬁtting using spectrum, the risk analysis using the convex Gaussian minimax
theorem, and the statistical inference using approximate message propagation methods.

1. Outline
1.1. Linear Regression. We consider a linear regression problem with 𝑝-dimensional covariates
𝑛 of a covariate 𝑋 ∈ R 𝑝 and a
and a parameter. Suppose that we observe i.i.d. 𝑛 pairs {(𝑋𝑖 , 𝑌𝑖 )}𝑖=1
𝑖target variable 𝑌𝑖 ∈ R generated from the following linear model with the true parameter 𝜃 0 ∈ R 𝑝 :𝑌𝑖 = 𝑋𝑖 , 𝜃 0  + 𝜉𝑖 , 𝑖 = 1, ..., 𝑛,where 𝜉𝑖 is a centered noise variable. Let Σ = E[𝑋𝑖 𝑋𝑖 ] be a covariance matrix of the covariate.
The goal of this problem is to estimate the parameter 𝜃 0 from the observations. Here, we consider
the high-dimensional setting, speciﬁcally, we consider 𝑝  𝑛 or where 𝑝 = ∞ regardless of 𝑛.
Also, we do not impose the assumption of sparsity on the true parameter 𝜃 0 . In this setting, thenotion of benign overﬁtting is actively studied.
We investigate whether benign overﬁtting-like phenomena occur in situations in which we relax
key assumptions in this foundational model. [NI22] considers the case where the covariates are
dependent in the sample direction and shows that similar benign overﬁtting occurs depending on
the strength of the dependence. [WI23] develops an informative prior distribution that allows
for a Bayesian estimator and its distribution approximation that is valid even in non-sparse highdimensions. [TI23] considers that the noise 𝜉𝑖 is dependent on the covariance and gives conditions
for the estimation error convergence, by utilizing the Gaussian minimax inequality.
1.2. Generalized Linear Regression. We next consider the generalized linear model (GLM):
for a pair (𝑋, 𝑌 ) of 𝑝-dimensional random features 𝑋 and random responses 𝑌 , we consider thefollowing modelE[𝑌 | 𝑋 = 𝑥] = 𝑔(𝑥  𝛽), ∀𝑥 ∈ R 𝑝 ,1(1)

where 𝑔 : R → R is an inverse link function that monotonically increases, and 𝛽 = (𝛽1 , . . . , 𝛽 𝑝 )  ∈𝑛
R 𝑝 is an unknown deterministic coeﬃcient vector. Suppose that we observe i.i.d. 𝑛 pairs {(𝑋𝑖 , 𝑌𝑖 )}𝑖=1
of a feature vector 𝑋𝑖 ∈ R 𝑝 and a target variable 𝑌𝑖 ∈ Y that follow the GLM (1), where Y is a
response space, such as R, R+ , {0, 1}, {0, 1, 2, . . .}, and so on.
We consider the proportional high-dimensional regime: we are particularly interested in the
proportional limit of the coeﬃcient dimension 𝑝 and sample size 𝑛:𝑛, 𝑝 → ∞ and 𝑝/𝑛 → ∃𝜅 ∈ (0, ∞).
In this regime, for the logistic regression as the special case, statistical inference on 𝛽 has been
actively studied without the sparsity of 𝛽.
[SUI23] develops a methodology for statistical inference for a broad class of GLMs, by deriving
the asymptotic normality of an estimator. This method is based on the analysis of a state evolution
equation by a vector approximate massage passiong (VAMP) and its application to statistical models.
1.3. Bandit with Linear Context. We study a bandit problem with 𝐾 arms associated with a linear
model for its rewards with 𝑝-dimensional context vectors. For each round 𝑡 ∈ [𝑇] := {1, 2, ..., 𝑇 }and arm 𝑖 ∈ [𝐾], we deﬁne a context 𝑋𝑡(𝑖) which is a 𝑝-dimensional zero-mean sub-Gaussian vector,
which is independent among rounds 𝑡. An agent chooses an arm 𝐼 (𝑡) ∈ [𝐾] based on 𝑋𝑡(𝑖) of all
the arms 𝑘 ∈ [𝐾], and then observes a reward that follows a linear model as shown in
𝑌 (𝐼 (𝑡)) = 𝑋 (𝐼 (𝑡)) , 𝜃 𝐼 (𝑡)  + 𝜉 (𝑡).
The unknown true parameters 𝜃 (𝑖) for each arm 𝑖 ∈ [𝐾] lie in a parameter space R 𝑝 , and the
independent sub-Gaussian noise 𝜉 (𝑡) with zero mean and variance 𝜎 2 > 0. We deﬁne 𝑖 ∗ (𝑡) :=argmax𝑖∈[𝐾] 𝑋𝑡(𝑖) , 𝜃 (𝑖)  as the (ex ante) optimal arm at round 𝑡.
Our goal is to design an algorithm that maximizes the total reward, which is equivalent tominimizing the following expected regret. [KI23] considers the high-dimensional setting 𝑝  𝑛 or
𝑝 = ∞ without the sparsity, then derive a novel explore-then-commit strategy to achieve minimize
the expected regret.
References[KI23] Junpei Komiyama and Masaaki Imaizumi. High-dimensional contextual bandit problem
without sparsity. Advances in Neural Information and Processing Systems, 2023.
[NI22] Shogo Nakakita and Masaaki Imaizumi. Benign overﬁtting in time series linear modelwith over-parameterization. arXiv preprint arXiv:2204.08369, 2022.
[SUI23] Kazuma Sawaya, Yoshimasa Uematsu, and Masaaki Imaizumi. Feasible adjustments ofstatistical inference in high-dimensional generalized linear models. 2023.
[TI23] Toshiki Tsuda and Masaaki Imaizumi. Benign overﬁtting of non-sparse high-dimensional
linear regression with correlated noise. arXiv preprint arXiv:2304.04037, 2023.
[WI23] Tomoya Wakayama and Masaaki Imaizumi. Bayesian analysis for over-parameterized
linear model without sparsity. arXiv preprint arXiv:2305.15754, 2023.
2

Statistical Challenges to Dimensionality in Astronomical Big Data
Tsutomu T. TAKEUCHI1,2, Kazuyoshi YATA3, Makoto AOSHIM3, Kohji YOSHIKAWA3, Kento EGASHIRA4,
Aki ISHII4, Suchetha COORAY5, Kouichiro NAKANISHI5, Kotaro KOHNO6, Ryusei R. KANO1, Yoh-Ichi
MOTOTAKE7, Taisei YAMAGATA1, Aina May SO1,8, Hai-Xia MA1, Shunya UCHIDA1, Wen E. SHI1,
Sena A. MATSUI1, and Kai T. KONO1
1. Nagoya University, Japan, 2. Institute of Statistical Mathematics, Japan, 3. University of Tsukuba, Japan, 4. Tokyo
University of Science, Japan, 5. National Astronomical Observatory of Japan, 6. The University of Tokyo, Japan, 7.
Hitotsubashi University, Japan, 8. Gakushuin University, Japan

1. Introduction: Galaxy Formation and Evolution with Big Data
Matter in the early Universe was almost uniform, and a slightly dense region grew by gravity, finally into a galaxy.
It was attempted to develop a theory to deal with the star formation and associated history of heavy element synthesis,
under an assumption that a galaxy has formed from a single, huge gas cloud. While the research in this direction
was once completed in the first half of 1980s, this was not the end of the studies of galaxy evolution. Cosmologicalresearch that has progressed in parallel has revealed that galaxies merge and grow. This indicates that the galaxy
evolution is a very complicated process that strongly depends on the density of the surrounding galaxies and the gas
density. In order to formulate the galaxy evolution, it is necessary to determine such a huge system of equations.
Though astrophysicists have constructed the governing equations from the physical laws from the first principle
before, such a method is not realistic anymore when the quantity space exceeds 10 dimensions. It is a high time torevolutionize the methodology for galaxy evolution studies.

2. Application of High-Dimensional Statistical Analysis to Astrophysical data
In astronomy, if we denote the dimension of data as d and the number of samples as n, we often meet a case withn << d. Traditionally, such a situation is regarded as ill-posed, and there was no choice but to throw away most of
the information in data dimension to let d < n. The data with n << d is referred to as high-dimensional low sample
size (HDLSS). To deal with HDLSS problems, a method called high-dimensional statistics has been developed
rapidly in the last decade. We first introduce the high-dimensional statistical analysis to the astronomical community.
We apply two representative methods in the high-dimensional statistical analysis methods, the noise-reduction
principal component analysis (NRPCA) and automatic sparse principal component analysis (A-SPCA), to aspectroscopic map of a nearby archetype starburst galaxy NGC 253 taken by the Atacama Large
Millimeter/Submillimeter Array (ALMA). The ALMA map is a typical HDLSS dataset. First, we analyzed theoriginal data including the Doppler shift due to the systemic rotation. The high-dimensional PCA could describe thespatial structure of the rotation precisely. We then applied to the Doppler-shift corrected data to analyze more subtle
spectral features. The NRPCA and R-SPCA could quantify the very complicated characteristics of the ALMA
spectra. Particularly, we could extract the information of the global outflow from the center of NGC 253. Thismethod can also be applied not only to spectroscopic survey data, but also any type of data with small sample sizeand large dimension. We are also trying to develop a method to analyze absorption line systems in the spectra of

distant radio quasars.

3. Galaxy Manifold
From 1970s to the mid-1980s, classical multivariate analysis methods such as the principal component analysis(PCA) were used to combine physical quantities of galaxies in a high-dimensional space. Various (logarithmic)
linear relations, so-called galactic scaling relations, have been discovered. Research to unify the scaling relations
and find the fundamental relationships has led to the concept of galaxy manifolds. However, the galaxy manifold
has once been almost forgotten because the classical PCA could treat only linear relations, and it remained a limited
concept, though they are still useful for exploring (log)linear relations of galaxies.
Recently, we discovered a galaxy manifold that expresses the basics of galactic evolution by the Fisher EM
algorithm. Because of its strongly nonlinear spatial structure, it could have never been found in previous studies
based on the classical PCA. To understand the manifold, a more sophisticated method beyond a mere classificationis needed. We focused on a method known as the manifold learning, one of the latest methods of data science thatis completely different from conventional methodologies.
We adopt the algorithm Isomap and UMAP (Uniform Manifold Approximation and Projection). Isomap defines
the neighboring points by using input-space distance and the distant points as a sequence of “short hops” between
neighboring points. Isomap tries to find shortest paths in a graph with edges connecting neighboring data points. By
construction, Isomap preserves the “surface density” of data points in the feature space. UMAP is based on
differential geometry and algebraic topology. The algorithm is founded on three assumptions: 1) the data are
uniformly distributed on a Riemannian manifold, 2) the Riemannian metric is locally constant (or can beapproximated as such), and 3) the manifold is locally connected. From these assumptions it is possible to model the
manifold with a fuzzy topological structure. Sine it defines the manifold so that the data points distribute ashomogeneously as possible, it does not preserve the surface density of data points. UMAP also preserves some
important structural properties, and it is more robust against noise than Isomap. Manifold learning algorithm can“unfold” a curved and/or rolled manifold in the feature space, and provide a local coordinate system on it. The
resulting manifolds with local coordinates from Isomap and UMAP are presented in Fig. 1. From Figure 1, weclearly see that the galaxy manifold is two-dimensional. We also stress that two different algorithms, Isomap and
UMAP yield similar two-dimensional manifolds. Since Isomap preserves the density of data point cloud, we observethat the manifold has a density structure, i.e., dense and sparse regions on the manifold.
The galaxy manifold obtained with Isomap preserve this information and reveal the speed of galaxy evolution at
various stages along the manifold. e.g., galaxies passes the green valley very fast. In contrast, the galaxy manifold
obtained with UMAP is imposed uniformity on the galaxy data, leading to a more robust and representativedescription of the observed galaxy properties e.g., galaxies evolve continuously in the feature space, without adiscontinuity or “jump” on their evolutionary tracks. Thus, the galaxy manifold provides a clue to the evolutionary
path of galaxies on the manifold. The SFR and stellar mass fields do not show the same evolutionary path. This
supports that the galaxy merger without star formation plays a significant role in the growth of stellar mass. Next
step is to fully parametrize the evolution equation of galaxies.

Predictive Density Estimation for Two Ordered Normal Means
Under α-Divergence Loss
Yuan-Tsung Chang (Mejiro University)
Nobuo Shinozaki (Keio University)
William, E. Strawderman (Rutgers University)

When the underlying loss metric is α-divergence, D(α), loss introduced by siszàr (1967), we
consider stochastic and Pitman closeness domination in predictive density estimation problems
when there are restrictions given on two means. The underlying distributions considered arenormal location-scale models, including the distribution of the observables, the distribution of
the variable whose density is to be predicted, and the estimated predictive density which will
be taken to be of the plug-in type. The scales may be known or unknown. The main contents
are as follows:
1. First, we introduce a general expression which derived by Chang and Strawderman (2014)
for the α-divergence loss as following:
If the true density function of Y is N (μ, σ 2 ) and the estimated predictive density of Y , is
N (μ̂, σ̂ 2 ) thena) for −1 < α < 1,(μ̂−μ)2
4
2
2 −A(σ 2 ,σ̂ 2 ) 2
Dα (N (ỹ|μ̂, σ̂ ), N (ỹ|μ, σ )) =
1
−
d(σ
,
σ̂
)e
,
1 − α222

whered(σ 2 , σ̂ 2 ) =σ (α−1)/2 τ,σ̂ (α+1)/2

A(σ 2 , σ̂ 2 ) =



1−α
2σ 2



1−

(1 − α)τ 2
2σ 2



> 0,

1
=
τ2



Further, d(σ 2 , σ̂ 2 ) < 1 and A(σ 2 , σ̂ 2 ) > 0.
b) (Reverse KL)
D+1 (N (ỹ|μ̂, σ̂ 2 ), N (ỹ|μ, σ 2 )) =

1
2

 2
σ̂

σ

− log
2









σ̂ 2
(μ̂ − μ)2
−
1
+
.
σ2
σ2

c) (KL)1
D−1 (N (ỹ|μ̂, σ̂ ), N (ỹ|μ, σ )) =
2
2

2

 2
σ

σ2
(μ̂ − μ)2
−
log
−
1
+
.
σ̂ 2
σ̂ 2
σ̂ 2

Also note that in each case, the {D(α)} loss is a concave monotone function of squared
error loss |μ̂ − μ|2 and is also a function of the variances. In this set-up and show that it isa concave monotone function of quadratic loss, and also of the variances (predicand, andplug-in).
2. Next, we demonstrate D(α) stochastic domination and Pitman closeness of certain plug-in
predictive densities over others for the entire class of metrics simultaneously when ”usual”
stochastic domination and Pitman closeness holds in the related problem of estimating
two ordered means with respect to quadratic loss(Oono, Shinozaki (2005), Chang, Oonoand Shinozaki(2012), Chang, Fukuda and Shinozaki(2017)).

11+α 1−α
+
.
2σ̂ 2
2σ 2

3. We also discuss improving the generalized Bayesian predictive densities suggested by Corcuera and Giummole (1999) under D(α) loss.
Based on the data Xij ∼ N (μi , σi2 ), i = 1, 2, j = 1, · · · , ni , we predict the density Ỹ ∼N (μi , σi2 ), i = 1, 2. We denote its density function by p(ỹ; μi , σi ), where μi and σi2 are
unknown.
When −1 ≤ α < 1, Corcuera and Giummole (1999) have established that the best invariant
predictive density of p(ỹ; μi , σi ) based solely on xi1 , · · · xini isp̂α (ỹ; x̄i , σ̃i ) ∝ 1 +y − x̄i1−α2ni + 1 − ασ̃i

2 −(2ni −1−α)/2(1−α),

where x̄i is the sample mean and σ̃i2 = ((ni −1)/ni )s2i is the sample variance. Corcuera andGiummole (1999) have also shown that p̂α (ỹ; x̄i , σ̃i ) is the generalized Bayesian predictive
density for the prior density f (μi , σi ) ∝ 1/σi , 0 < σi < ∞. It is to be noted that p̂α (ỹ; x̄i , σ̃i )is not a normal distribution, although the plug-in density N (x̄i , s2i ) is the generalized Bayes
rule when α = 1.
We consider the following two cases separately where order restrictions on μi and/or σi2
are present,i) Case when μ1 ≤ μ2 .
ii) Case when μ1 ≤ μ2 and σ12 ≤ σ22 .
Examples of D(α) stochastic (Pitman closeness) domination presented relate to the problem of estimating the predictive density of the variable with the restrictions on two normalmeans.
keywords: Predictive density, α-divergence, stochastic dominance, ordered normalmeans, Pitman closeness criterion

References[1] Aitchison, J., Goodness of prediction ﬁt, Biometrika, 62 (1975) 545-554.
[2] Chang, Y.-T., Oono, Y., Shinozaki, N., Improved estimators for the common mean andordered means of two normal distributions with ordered variances. Journal of Statistical
Planning and Inference, 142 (2012) 2619-2628.
[3] Chang, Y.-T., Strawderman, W. E., Stochastic domination in predictive density estimationfor ordered normal means under α-divergence loss. Journal of Multivariate Analysis, 128(2014) 1-9.
[4] Chang, Y.-T., Fukuda, K., Shinozaki, N., Estimation of two ordered normal means when acovariance matrix is known. Statistics, 5 (2017) 1095-1104.
[5] Chang, Y.-T., Shinozaki, N, Strawderman, W. E., Pitman closeness domination in predictive density estimation for two-ordered normal means under α-divergence loss. Japanese
Journal of Statistics and Data Science, (2020) 3:1-21.
[6] Csiszàr, I., Information-type measures of diﬀerence of probability distributions and indirect
observations, Studia Sci. Math. Hungar., 2 (1967) 299-318.
[7] Fourdrinier, D., Marchand, E., Righi, A., Strawderman, W. E., On Improved predictive
density with parametric constraints. Electronic Journal of Statistics, 5 (2011) 172-191.
[8] Oono, Y., Shinozaki, N., Estimation of two order restricted normal means with unknown
and possibly unequal variances. Journal of Statistical Planning and Inference, 131 (2005)(2), 349-363.

2

On the eﬃciency-loss free ordering-robustness of
product-PCA
Hung Hung
1 Institute of Health Data Analytics and Statistics, National Taiwan University, Taiwan

Abstract
This article studies the robustness of the eigenvalue ordering, an important issue
when estimating the leading eigen-subspace by principal component analysis (PCA).
In Yata and Aoshima (2010), cross-data-matrix PCA (CDM-PCA) was proposed andshown to have smaller bias than PCA in estimating eigenvalues. While CDM-PCA
has the potential to achieve better estimation of the leading eigen-subspace than theusual PCA, its robustness is not well recognized. In this article, we ﬁrst develop a
more stable variant of CDM-PCA, which we call product-PCA (PPCA), that provides
a more convenient formulation for theoretical investigation. Secondly, we prove that,in the presence of outliers, PPCA is more robust than PCA in maintaining the correct
ordering of leading eigenvalues. The robustness gain in PPCA comes from the random
data partition, and it does not rely on a data down-weighting scheme as most robust
statistical methods do. This enables us to establish the surprising ﬁnding that, when
there are no outliers, PPCA and PCA share the same asymptotic distribution. Thatis, the robustness gain of PPCA in estimating the leading eigen-subspace has noeﬃciency loss in comparison with PCA. Simulation studies and a face data example
are presented to show the merits of PPCA. In conclusion, PPCA has a good potential
to replace the role of the usual PCA in real applications whether outliers are present
or not.
Key words: cross-data-matrix PCA; dimension reduction; eﬃciency loss; orderingof eigenvalues; random partition; robustness.

1/HDUQLQJ2UGLQDOLW\LQ+LJK'LPHQVLRQDO'DWD-HRQJ\RXQ$KQ
'HSDUWPHQWRI,QGXVWULDODQG6\VWHPV(QJLQHHULQJ.$,67
1XPHURXVUHDOZRUOGDSSOLFDWLRQVLQYROYHQDWXUDOO\RUGLQDORXWFRPHVVXFKDVFDQFHUVWDJHVRUWXPRUJUDGHV'HVSLWHWKHUHFHQWVXUJHLQKLJKGLPHQVLRQDOVWDWLVWLFDOPHWKRGRORJLHVKLJKGLPHQVLRQDOOHDUQLQJZLWKRUGLQDORXWFRPHVKDVEHHQODUJHO\RYHUORRNHGLQWKH+'/66OLWHUDWXUH,QWKLVWDON,ZLOOLQWURGXFHUHFHQWSURMHFWVRQRUGLQDOLW\LQKLJKGLPHQVLRQDOGDWD7KHILUVWWKUHHWRSLFVFRQFHUQVXSHUYLVHG
OHDUQLQJDLPHGDWSUHGLFWLQJRUGLQDOODEHOV$OOWKUHHRUGLQDOPHWKRGVDVVXPHVSDUVLW\DQGHTXDOFRYDULDQFHSRSXODWLRQVWUXFWXUHOHDGLQJXVWRWHUPWKHP RUGLQDOVSDUVHKLJKGLPHQVLRQDO/'$ 7KH\RSHUDWHRQWKHSULQFLSOHWKDWDFODVVLILFDWLRQUXOHSULPDULO\UHOLDQWRQYDULDEOHVWKDWDUHPRQRWRQLFDOO\DVVRFLDWHGZLWKWKH
RXWFRPHVKRXOGEHSUHIHUDEOH7KH\DOOUHVXOWLQDORZGLPHQVLRQDOGLVFULPLQDQWVXEVSDFHZKHUHFODVVHVDUHVHTXHQWLDOO\DOLJQHG7KHILUVW):2&PHWKRGZHLJKWVIHDWXUHVEDVHGRQWKHLUUDQNFRUUHODWLRQVZLWKFODVVODEHOVLQWHJUDWLQJWKHVHZHLJKWVLQWRWKH/'$IUDPHZRUN7KHVHFRQG62%/PHWKRGFRPELQHVVSDUVLW\DQG
RUGLQDOLW\UHJXODUL]DWLRQVLQDKLJKGLPHQVLRQDOJHQHUDOL]HGHLJHQYDOXHSUREOHP7KH
WKLUG62'$DSSURDFKDSSOLHGUHJXODUL]DWLRQWRRSWLPDOVFRUHVZLWKLQWKHVSDUVH/'$
IUDPHZRUN,QDGGLWLRQLQVFHQDULRVZKHUHDQRUGLQDORXWFRPHLVXQREVHUYHGRQHPD\VHDUFKIRUDQRUGLQDOVLJQDOLQWKHGDWD7KLVOHDGVXVWRGHYHORS PRQRWRQH
FOXVWHULQJ WKDWLVGHVLJQHGWRLGHQWLI\VXEJURXSVLQWHUSUHWDEOHLQDQRUGLQDO
PDQQHUNormal-reference test for high-dimensional covariance matrices
Jin-Ting Zhang
National University of Singapore

Abstract
In the past decade, much attention has been paid for testing the equality of high-dimensional
covariance matrices. Several test statistics have been proposed for this purpose. Some of them
imposedstrong assumptions, aiming to yield the asymptotic normality of the associated test

statistics. In practice, however, these assumptions are often challenging to verify, resulting in size
control issues when the required assumptions are not met. To address this challenge, in this talk,
we investigate anormal-reference test which can effectively control the size. In the normal-

reference test, the null distribution of a test statistic is approximated with that of a chi-square-type
mixture which is obtained from the test statistic under the null hypothesis, assuming normality of thedata samples. To accurately approximate the distribution of the chi-square-type mixture, we employ
a three-cumulant matched ǘ2-approximation with the approximation parameters being consistently
estimated from the data. Two simulation studies demonstrate that in terms of size control, theproposed normal- reference test performs well across a range of scenarios and it outperforms several
existing competitors. A real data example illustrates the proposed normal-reference test.

KEY WORDS: ǘ2-type mixtures; high-dimensional data; three-cumulant matched ǘ2-approximation.

   # !    "    % "((#!$!""#! %"" (4 )!"#,!/)"6 #$"$"$#(#/#5% "$3"!!, #""#!(" #!"!  %  "! # !  #"$ ! , #"" ""# *#  3"($)!#!!!""(!#"&*!#"#,#!" "" !##" "-/#"#,# !#!""#2"""($"))!$#"#"!"( !""(!!"#""*!*##"", #$!*!2 ! "!(!-$""##,# !$ "#$"$" ,  , ! "! # # )("  #  !)!#!+#!!""!"("2 ! "#*!#"""!(!-#"#"#*#!##, ""#!(#(!""( $"#)!#!+#"#!!!""145#" #!"(!#")!)!"##!)#045#")!"" )!"#!(#(!2"*##"/# ! "#"#""#, !)# !!#!$#"#2  " !"" # !   #  $ !(!-$ !#! *# "3#!$!*!, $ !"$!($##!$)"2" $/*"!# !#$ """"$"()!"(!#")(#!"(!#"!)!("!!"2
4 ""#*!*#! /+!( 52PL SS_PsAJh1 aJSGAL: 6_PJ LPL@GP:@*PL*o1
LPL@aJPPh> .Aah_A"lhAPLa oA  GL:1oAL@huS1 JPLh1*_GP G:P_Ah>J
a>P:P LEEAh

q2 +QMbB/2` i?2 T`Q#H2K Q7 bKTHBM; 7`QK  :B##b /Bbi`B#miBQM π(/x) ∝ 2tT(−U (x))/x
QM (_d , B(_d ))- r?2`2 U : _d → [0, ∞) Bb  MQM@M2;iBp2 TQi2MiBH 7mM+iBQMX PM2 Q7 i?2
2ti2MbBp2Hv mb2/ ivT2b Q7 H;Q`Bi?Kb 7Q` i?2 bKTHBM; Bb i?2 GM;2pBM ivT2 KQiBpi2/ #v
i?2 GM;2pBM /vMKB+b- i?2 bQHmiBQM Q7 i?2 7QHHQrBM; d@/BK2MbBQMH biQ+?biB+ /Bz2`2MiBH
2[miBQM Ua.1V,√
UyXRV/Xt = −∇U (Xt ) /t + 2/Bt , X0 = ξ,r?2`2 {Bt }t≥0 Bb  d@/BK2MbBQMH "`QrMBM KQiBQM M/ ξ Bb  d@/BK2MbBQMH `M/QK
p2+iQ` rBi? |ξ| < ∞ HKQbi bm`2HvX aBM+2 i?2 k@qbb2`bi2BM Q` iQiH p`BiBQM /BbiM+2#2ir22M π M/ i?2 Hr Q7 Xt Bb +QMp2`;2Mi mM/2` KBH/ +QM/BiBQMb- r2 2tT2+i i?i i?2
Hrb Q7 GM;2pBM@ivT2 H;Q`Bi?Kb BMbTB`2/ #v Xt b?QmH/ +QMp2`;2 iQ πX >Qr2p2`- KQbi Q7i?2 i?2Q`2iB+H ;m`Mi22b 7Q` bm+? H;Q`Bi?Kb `2 #b2/ QM i?2 +QMp2tBiv Q7 U - i?2 irB+2+QMiBMmQmb /Bz2`2MiB#BHBiv Q7 U - Q` i?2 GBTb+?Bix +QMiBMmBiv Q7 i?2 ;`/B2Mi ∇U - r?B+? /Q
MQi ?QH/ BM bQK2 KQ/2HHBM; BM biiBbiB+b M/ K+?BM2 H2`MBM;X h?2 KBM BMi2`2bi Q7 i?Bb
bim/v Bb T`QTQbH Q7  GM;2pBM@ivT2 H;Q`Bi?K r?Qb2 +QMp2`;2M+2 +M #2 ;Bp2M mM/2`
KBMBKH bbmKTiBQMbX
hQ b22 r?i /B{+mHiB2b r2 M22/ iQ /2H rBi?- r2 `2pB2r  ivTB+H MHvbBb (e) #b2/ QM
i?2 bKQQi?M2bb Q7 U - i?i Bb- i?2 irB+2 +QMiBMmQmb /Bz2`2MiB#BHBiv Q7 U M/ i?2 GBTb+?Bix+QMiBMmBiv Q7 ∇U X 6B`biHv- i?2 irB+2 +QMiBMmQmb /Bz2`2MiB#BHBiv bBKTHB}2b /Bb+mbbBQMb Q`
THvb bB;MB}+Mi `QH2b BM bim/B2b Q7 7mM+iBQMH BM2[mHBiB2b bm+? b SQBM+`û BM2[mHBiB2b
M/ HQ;`Bi?KB+ aQ#QH2p BM2[mHBiB2b (2X;X- R- k)X aBM+2 i?2 7mM+iBQMH BM2[mHBiB2b 7Q` π
`2 2bb2MiBH BM MHvbBb Q7 GM;2pBM H;Q`Bi?Kb- i?2 bbmKTiBQM i?i U Bb Q7 +Hbb C 2
7`2[m2MiHv TT2`b BM T`2pBQmb bim/B2bX AM i?2 b2+QM/ TH+2- i?2 GBTb+?Bix +QMiBMmBiv+QK#BM2/ rBi? r2F +QM/BiBQMb 2Mbm`2b i?2 `2T`2b2MiiBQM Q7 i?2 HBF2HB?QQ/ `iBQ #2ir22M{Xt } M/ {Yt }- r?B+? Bb +`BiB+H r?2M r2 #QmM/ i?2 EmHH#+FĜG2B#H2` /Bp2`;2M+2X GBTib2`
M/ a?B`v2p (9) 2t?B#Bi Km+? r2F2` +QM/BiBQMb i?M LQpBFQpǶb Q` ExKFBǶb +QM/BiBQM
7Q` i?2 2tTHB+Bi `2T`2b2MiiBQM B7 UyXRV ?b i?2 mMB[m2 bi`QM; bQHmiBQMX aBM+2 i?2 GBTb+?Bix+QMiBMmBiv Q7 ∇U Bb bm{+B2Mi 7Q` i?2 2tBbi2M+2 M/ i?2 mMB[m2M2bb Q7 i?2 bi`QM; bQHmiBQM
Q7 UyXRV- i?2 7`K2rQ`F Q7 GBTib2` M/ a?B`v2p (9) Bb TTHB+#H2X
Pm` TT`Q+?2b iQ Qp2`+QK2 i?2 MQM@bKQQi?M2bb Q7 U `2 KQHHB}+iBQM-  +HbbB+H
TT`Q+? iQ /2HBM; rBi? MQM@bKQQi?M2bb BM /Bz2`2MiBH 2[miBQMb- M/ i?2 ǵKBbmb2Ƕ Q7
KQ/mHB Q7 +QMiBMmBiv 7Q` TQbbB#Hv /Bb+QMiBMmQmb 7mM+iBQMbX q2 +QMbB/2` i?2 +QMpQHmiBQM
Ūr := U ∗ ρr QM U rBi?  r2F ;`/B2Mi- M/ bQK2 bm{+B2MiHv bKQQi? MQM@M2;iBp2
7mM+iBQM ρr rBi? +QKT+i bmTTQ`i BM  #HH Q7 +2Mi`2 y M/ `/Bmb r ∈ (0, 1]X q2 +M
H2i Ūr #2 Q7 +Hbb C 2 M/ Q#iBM #QmM/b 7Q` i?2 +QMbiMi Q7 SQBM+`û BM2[mHBiB2b 7Q`π̄ r (/x) ∝ 2tT(−Ūr (x))/x- r?B+? bm{+2 iQ b?Qr i?2 +QMp2`;2M+2 Q7 i?2 Hr Q7 i?2 KQHHB}2/
EQK# AMbiBimi2 7Q` a+B2M+2- lMBp2`bBiv Q7 hQFvQ- j@3@R EQK#- J2;m`Q@Fm- hQFvQ
R8j@3Nyk- CTM
1@KBH //`2bb, MFFBi!;X2++Xm@iQFvQX+XDTX
h?2 mi?Q` rb bmTTQ`i2/ #v CaSa EE1L>A :`Mi LmK#2` CSkREkyjR3 M/ Cah *_1ah :`Mi
LmK#2`b CSJC*_kR.k M/ CSJC*_kRR8X
R

ka LEEAh

/vMKB+b {X̄tr } /2}M2/ #v i?2 a.1√ /X̄tr = −∇Ūr X̄tr /t + 2/Bt , X̄0r = ξ

iQ i?2 +Q``2bTQM/BM; :B##b /Bbi`B#miBQM π̄ r BM k@qbb2`bi2BM /BbiM+2 QrBM; iQ "F`v
2i HX (R)- GBm (8)- M/ G2?2+ (j)X aBM+2 i?2 +QMpQHmiBQM ∇Ūr Bb GBTb+?Bix +QMiBMmQmb
B7 i?2 KQ/mHmb Q7 +QMiBMmBiv Q7  `2T`2b2MiiBp2 ∇U Bb }MBi2 Ui?2 +QMp2`;2M+2 iQ x2`Q
Bb mMM2+2bb`vV-  +QM+Bb2 `2T`2b2MiiBQM Q7 i?2 HBF2HB?QQ/ `iBQb #2ir22M i?2 KQHHB}2//vMKB+b {X̄tr } M/ {Yt } Bb pBH#H2- M/ r2 +M 2pHmi2 i?2 EmHH#+FĜG2B#H2` /Bp2`;2M+2
mM/2` r2F bbmKTiBQMbX
b Qm` MHvbBb `2HB2b QM KQHHB}+iBQM- i?2 #BbĜp`BM+2 /2+QKTQbBiBQM BM 2biBKiBQM Q7∇Ūr `i?2` i?M ∇U Bb +`m+BHX h?Bb /2+QKTQbBiBQM 2M#H2b mb iQ T`QTQb2 M2r H;Q`Bi?Kb
7Q` U rBi?Qmi +QMiBMmQmb /Bz2`2MiB#BHBivX *QM+`2i2Hv bT2FBM;- r2 T`QTQb2  M2r H@;Q`Bi?K MK2/ i?2 bT?2`B+HHv bKQQi?2/ GM;2pBM JQMi2 *`HQ Uaa@GJ*V H;Q`Bi?Kr?Qb2 2``Q`b +M #2 `#Bi``BHv bKHH mM/2` i?2 /BbbBTiBpBiv Q7 U M/ i?2 #QmM/2/M2bb Q7
i?2 KQ/mHmb Q7 +QMiBMmBiv Q7 r2F ;`/B2MibX AM //BiBQM- r2 `;m2 x2`Qi?@Q`/2` p2`bBQMb
Q7 i?2b2 H;Q`Bi?Kb r?B+? `2 Mim`HHv Q#iBM2/ pB BMi2;`iBQM #v T`ibX
_272`2M+2b(R) "F`v- .X- "`i?2- 6X- *iiBmt- SX- M/ :mBHHBM- X Ukyy3VX  bBKTH2 T`QQ7 Q7 i?2
SQBM+`û BM2[mHBiv 7Q`  H`;2 +Hbb Q7 T`Q##BHBiv K2bm`2bX 1H2+i`QMB+ *QKKmMB+iBQMb
BM S`Q##BHBiv- Rj,eyĜeeX
(k) *iiBmt- SX- :mBHHBM- X- M/ qm- GX@JX UkyRyVX  MQi2 QM hH;`M/ƪ
b i`MbTQ`iiBQM
BM2[mHBiv M/ HQ;`Bi?KB+ aQ#QH2p BM2[mHBivX S`Q##BHBiv h?2Q`v M/ _2Hi2/ 6B2H/bR93,k38Ĝjy9X(j) G2?2+- CX UkykRVX h?2 GM;2pBM JQMi2 *`HQ H;Q`Bi?K BM i?2 MQM@bKQQi? HQ;@+QM+p2+b2X hQ TT2` BM i?2 MMHb Q7 TTHB2/ S`Q##BHBivX(9) GBTib2`- _X aX M/ a?B`v2p- X LX UkyyRVX aiiBbiB+b Q7 _M/QK S`Q+2bb2b, AX :2M2`H
i?2Q`vX aT`BM;2`- kM/ 2/BiBQMX(8) GBm- uX UkykyVX h?2 TQBM+`û BM2[mHBiv M/ [m/`iB+ i`MbTQ`iiBQM@p`BM+2 BM@
2[mHBiB2bX 1H2+i`QMB+ CQm`MH Q7 S`Q##BHBiv- k8URV,RĜReX
(e) _;BMbFv- JX- _F?HBM- X- M/ h2H;`bFv- JX UkyRdVX LQM@+QMp2t H2`MBM; pB
biQ+?biB+ ;`/B2Mi GM;2pBM /vMKB+b,  MQMbvKTiQiB+ MHvbBbX AM S`Q+22/BM;b Q7
i?2 kyRd *QM72`2M+2 QM G2`MBM; h?2Q`v- T;2b Red9ĜRdyjX

Two step estimations via the Dantzig selector for ergodic
time series models
Kou Fujimori1 and Koji Tsukuda2
1

Faculty of Economics and Law, Shinshu University.
2
Faculty of Mathematics, Kyushu University.

Let us consider the following time series model.
Xt = S(α φ(Xt−1 ), β  Zt−1 ) + ut ,2

2

E[ut |Ft−1 ] = σ (Xt−1 ; h),where Xt−1 = (Xt−1 , . . . , Xt−d ), without loss of generality. Let θ0 = (α0 , β0 ) be thetrue value of θ = (α , β  ) , Θ = Θα × Θβ ⊂ Rp+q a parameter space for θ and H a metricspace equipped with a metric dH . Put T10 := {j : α0j = 0}, T20 := {j + p1 : β0j = 0} andT0 = T10 ∪ T20 . We observe (X1 , Z1 ), . . . , (Xn , Zn ). Our aim is to estimate θ = (α , β  ) .
Hereafter, we ﬁx an initial value (X0 , . . . , X1−d ) = (x0 , . . . , x1−d ) and put p = p1 + p2 ,s = s1 + s2 , where s1 and s2 are the numbers of elements in T10 and T20 , respectively.
(1)
We ﬁrst construct the estimator θ̂n for θ by the following Dantzig selector type estimator:θ̂n(1) := arg min θ 1 , Cn := {θ ∈ Rp : ψn(1) (θ) ∞ ≤ λn },
θ∈Cn

where
1 ∂
S(α φ(Xt−1 ), β  Zt−1 ){Xt − S(α φ(Xt−1 ), β  Zt−1 )}n∂θnψn(1) (θ) =t=1

and λn is a tuning parameter. Moreover, we deﬁne the following estimator T̂n for T0 :T̂n := {j : |θ̂nj | > τn },where τn is a threshold. For the second step, we construct a consistent estimator for h, by(1)(2)
using θ̂n . Finally, using T̂n and ĥn , we consider the estimator θ̂n for θ as a solution to
the following equation:
ΨnT̂n (θT̂n , ĥn ) = 0, θ̃nT̂ c = 0,n1whereΨnT (θT , ĥn ) =∂n1  ∂θT S(αT1 φ(Xt−1 )T1 , βT2 Zt−1T2 )nσ 2 (Xt−1 ; ĥn )t=1·{Xt − S(αT1 φ(Xt−1 )T1 , β  Zt−1T2 )}for every T = T1 ∪ T2 .
(1)
In this talk, we establish the rate of convergence of θ̂n , and the asymptotic normality of θ̃nT̂n under some regularity conditions under high-dimensional and sparse settings.
Moreover, we discuss the integer-valued autoregressive models as an example.

Acknowledgements
This work was supported by JSPS KAKENHI Grant Numbers 21K13271(K.F.), 21K13836(K.T.).

2

Innovation algorithm of fractionally integrated (I(d)) process
and applications on the estimation of parameters
Junichi Hirukawa and Kou Fujimori
Niigata University and Shinshu University
ABSTRACTThe long memory phenomena frequently occur in the empirical studies of various ﬁelds. The fractionally integrated
process is the one of the suitable candidate which appropriately represents the long memory property. There are tworecursive algorithms for determining the one-step predictors of time series, that is, the Durbin-Levinson algorithm andthe innovation algorithm. The Durbin-Levinson algorithm for the fractionally integrated process is well-known andwidely used, which naturally derives the Cholesky factorization of the inverse matrix of the covariance matrix of theprocess. In this paper, we derive the innovation algorithm for the fractionally integrated process. The result is also
applied to the derivations of the Cholesky factorization of the covariance matrix and the Gaussian likelihood of theprocess in the explicit forms. Moreover, the asymptotic theory of Gaussian maximum likelihood estimator (GMLE) isderived in terms of the innovation algorithm.

1 Introduction
An ARMA (p,q) process {xt } is often called a short memory process since the covariance between xt and xt+ j decreases
rapidly as j → ∞. However, the long memory phenomena frequently occur in the empirical studies of various ﬁelds(see e.g., Hurst (1951)). In this paper, we consider one of the long memory process so-called the fractionally integrated(I(d)) process deﬁned by(1)(1 − L)d zt = εt ,(t = 1, . . . , n) ,

i.i.d.
where d ∈ (−1/2, 1/2), (d  0), L is the lag operator and {εt } ∼ 0, σ2 . Using the expansion of the lag operator(2)

Δ(L) = (1 − L)d =∞∞1  Γ( j − d) j 
ϕ j L j,L =Γ(−d) j=0 Γ( j + 1)
j=0

this can be rewritten as
(3)

εt =

∞


ϕ j L j zt =

j=0

∞


ϕ j zt− j .

j=0Then, {zt } is a stationary long memory process generated by
zt = (1 − L)−d εt∞∞
1  Γ( j + d)=
ψ j εt− j ,εt− j =
Γ(d) j=0 Γ( j + 1)
j=0
where the coeﬃcients satisfy ψ j = O jd−1 , so that the degree of decreasing is quite slow as j → ∞.

1

2 Main result
In this section, we provide the main results.

2.1 The Gaussian MLE for I(d) process
In this section, we impose the Gaussian assumption on I(d) process. Then, we have the Gaussian log-likelihood ofI (d) process for θ = d, σ2nn1  u j−1 (d)2n1.
l (θ) = l d, σ2 = − log {2π} −log v j−1 (θ) −2
2 j=1
2 j=1 v j−1 (θ)(4)

Now, we have the following main results. First, we describe the consistency of GMLE.
 = d,

Theorem 1. Let {zt } is the Gaussian I (d) process deﬁned in (1) with d ∈ (−1/2, 1/2), (d  0). And let θ
σ2 is


 has
the Gaussian MLE (GMLE) of θ = d, σ2 which maximizes the Gaussian log-likelihood (4). Then, the GMLE θ
consistency, that is,
P
→
θ0 ,
θ



where θ0 = d0 , σ20 is the true value of θ.
Next, we have the following asymptotic normality.
 = d,
σ2Theorem 2. Let {zt } is the Gaussian I (d) process deﬁned in (1) with d ∈ (−1/2, 1/2), (d  0). And let θ
is the Gaussian MLE (GMLE) of θ = d, σ2 which maximizes the Gaussian log-likelihood (4). Then, the GMLE θsatisﬁes the asymptotic normality, that is, L√  − θ0 →n θ
N (0, A) ,
where θ0 = d0 , σ20 is the true value of θ.

References
Akaike, H. (1969). Power Spectrum Estimation through Autoregressive Model Fitting. Ann. Inst. Stat. Math. 21,
407–419.
Berk, K. N. (1974). Consistent Autoregressive Spectral Estimates. Ann. Statist. 2, 489–502.
Brockwell, P. J. and Davis, R. A. (1991). Time Series: Theory and Methods. New York: Springer.
Hosking, J. R. M. (1981). Fractional Diﬀerencing. Biometrika 1, 165–176.
Hurst, H. E. (1951). Long-Term Storage Capacity of Reservoirs. Trans. Amer. Soc. Civil Eng. 116, 770–799.
Li, W. K. and McLeod, A. I. (1986). Fractional Time Series Modelling. Biometrika 73, 217-221.
Slater, L J. (1966). Generalized Hypergeometric Functions. New York: Cambridge University Press.
Yajima, Y. (1985). On Estimation of LongʖMemory Time Series Models. Australian Journal of Statistics 27, 303–320.
2

Scaling Limits of Markov Chains/Processes in Monte
Carlo Methods
Kengo Kamatani (Institute of Statistical Mathematics, JST CREST)

In this presentation, we will explore the recent results of scaling limit of piecewise deterministic
Markov processes for anisotropic targets. Suppose we wish to sample from
Πpdxq “ expp´Hpxqqdx
where H : Rd Ñ R is a continuously diﬀerentiable function. For the Bayesian context, thisprobability distribution is the posterior distribution of interest. If we have an i.i.d. sample fromΠ, we can approximate Π-integral of any function f pxq by the law of large numbers. In mostof the cases, direct i.i.d. sampling is impossible or computationally very expensive. For thesecases, the Markov chain Monte Carlo method is useful which originated with the classic paper byMetropolis et al. (1953) almost 70 years ago. The Markov chain Monte Carlo method is designed
to construct an ergodic Markov kernel P which is Π-invariant. If a Markov chain X1 , X2 , . . . isgeneraetd from the Markov kernel P then the law of large numbers is satisﬁed. The Markov
chain Monte Carlo is now a gold standard for Bayesian inference.
Recently, its continuous process version, the Markov process Monte Carlo method is ofsubstantial interest for Monte Carlo analysis. Known Markov process Monte Carlo methods relyon an auxiliary variable trick which uses an auxiliary variable v with a probability density ν
on Ξ and considers the joint probability distribution μ :“ Πpdxq b νpdvq as an extended target
distribution on Z “ Rd ˆ Ξ. The original target distribution is a marginal distribution of theextended target distribution. Since Brownian motion does not have an absolutely continuous
path, we can not simulate processes driven by Brownian motion exactly. For our Monte Carloanalysis, exact sampling is necessary. Therefore, the Markov processes of interest should nothave a Brownian part. Known processes consist of a deterministic part and a pure jump part.
These processes are known as the piecewise deterministic Markov processes.
Here we follow Azaïs et al. (2014) for the expression of the piecewise deterministic Markov
processes. The processes are constructed by characteristics pφ, λk , Qk : l “ 1, . . . , Kq. The ﬂow
φ : Z ˆ R Ñ Z is continuous, φp¨, tq is a homeomorphism for each t P R and φpφp¨, sq, tq “
φp¨, s ` tq. For each k “ 1, . . . , K, the jump rate λk : Z Ñ R` determines the jump time of pureştjump processes, and Qk is a Markov kernel on Z. Let Λk pz, tq “ 0 λk pφpz, sqqds.
The Markov process is deﬁned by the following way. Suppose zp0q “ pxp0q, tp0qq P Z. Let
T1 , . . . , TK be independent processes with PpTk ě tq “ expp´Λk pz, tqq. Let T˚ “ mink“1,...,K Tk .
If Tk “ T˚ , then Z is generated from Qk pφpz, T˚ q, ¨q and set"
φpzp0q, tq for t ă T˚
Xptq “Z
for t “ T˚ .
After T˚ , the process evolves in the same way with starting value Z. There are several choicesof characteristics. Two popular piecewise deterministic Markov processes use the same ﬂow φdeﬁned by x1 ptq “ vptq and v 1 ptq “ 0. The Zig-Zag sampler proposed by Bierkens et al. (2019)1

uses d Markov kernels Q1 , . . . , Qd with d jump rates λ1 , . . . , λd . For each i “ 1, . . . , d, the Markov
kernel is a deterministic kernel Qi deﬁned by a map px, vq ÞÑ px, Fi pvqq where Fi is an operator
that ﬂips the i-th coordinate of x. The jump rate is deﬁned by λi ppx, vqq “ maxt0, Bi Hpxqvi u.
The bouncy particle sampler proposed by Peters and de With (2012), Bouchard-Côtéet al. (2018) uses two Markov kernels Qbounce and Qref with corresponding jump rates λbounceand λref . The kernel Qbounce is a deterministic kernel deﬁned by a map px, vq ÞÑ px, κpx, vqq:κpx, vq “ v ´ 2x∇Hpxq, vy∇Hpxq}∇Hpxq}2and λbounce px, vq “ maxt0, x∇Hpxq, vyu. The jump rate λref is a positive constant, and Qref is a μinvariant Markov kernel. For our analysis, for simplicity, we assume Qref ppx, vq, dpy, wqq “ νpdwq.
We have several critical ﬁndings. For the Zig-Zag algorithm, its performance is intricately
linked to the orientation of the target’s anisotropy; speciﬁc alignments with the algorithm’s
operational axes lead to enhanced eﬃciency, while others can hinder its eﬀectiveness. The BPS
algorithm, on the other hand, exhibits a deterministic dynamical behaviour in its limiting form
with a better rate of convergence.
This is joint work with Joris Bierkens (TU Delft) and Gareth O. Roberts (Warwick). See ourpaper on arxiv https://arxiv.org/abs/2305.00694 for the detail.

References
Christophe Andrieu, Alain Durmus, Nikolas Nüsken, and Julien Roussel. Hypocoercivity ofpiecewise deterministic markov process-monte carlo. The Annals of Applied Probability, 31, 102021. ISSN 1050-5164. doi: 10.1214/20-AAP1653.
Romain Azaïs, Jean-Baptiste Bardet, Alexandre Génadot, Nathalie Krell, and Pierre-André Zitt.
Piecewise deterministic Markov process—recent results. In Journées MAS 2012, volume 44 of
ESAIM Proc., pages 276–290. EDP Sci., Les Ulis, 2014. doi: 10.1051/proc/201444017.
Joris Bierkens, Paul Fearnhead, and Gareth Roberts. The zig-zag process and super-eﬃcient
sampling for Bayesian analysis of big data. Ann. Statist., 47(3):1288–1320, 2019. ISSN 00905364. doi: 10.1214/18-AOS1715.
Joris Bierkens, Kengo Kamatani, and Gareth O. Roberts. High-dimensional scaling limits ofpiecewise deterministic sampling algorithms. The Annals of Applied Probability, 32, 10 2022.
ISSN 1050-5164. doi: 10.1214/21-AAP1762.
Alexandre Bouchard-Côté, Sebastian J. Vollmer, and Arnaud Doucet. The bouncy particle sampler: A nonreversible rejection-free markov chain monte carlo method. Journal of the American
Statistical Association, 113(522):855–867, 2018. doi: 10.1080/01621459.2017.1294075.
George Deligiannidis, Daniel Paulin, Alexandre Bouchard-Côté, and Arnaud Doucet. Randomized hamiltonian monte carlo as scaling limit of the bouncy particle sampler and dimensionfree convergence rates. The Annals of Applied Probability, 31, 12 2021. ISSN 1050-5164. doi:
10.1214/20-AAP1659.
N. Metropolis, W. A. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equations ofstate calculations by fast computing machines. Journal of Chemical Physics, 21:1087–1092,
1953.
Elias AJF Peters and G. de With. Rejection-free monte carlo sampling for general potentials.
Physical Review E, 85(2):026703, 2012.

On a general linear hypothesis testing problem for latent factor
models in high dimensions
Takahiro Nishiyamaa and Masashi HyodobaDepartment of Business Administration, Senshu Universityb
Faculty of Economics, Kanagawa University

Let xgi “ pxgi1 , . . . , xgip qJ „ Fg be iid p-dimensional random vectors collected fromthe ith subject in the gth population, where Fg denotes the distribution function forgth population, i P t1, . . . , ng u, g P t1, . . . , ku. A factor model assumes that for eachg P t1, . . . , ku, the observable vector xgi is decomposable into a latent factor and an
idiosyncratic component as follows:
xgi “ μg ` Fg zgi ` Ψ1{2g gi ,

(1)

where μg P Rp is a deterministic intercept vector, zgi “ pzgi1 , . . . , zgidg qJ is a dg -dimensional
latent factor vector, and gi “ pgi1 , . . . , gip qJ is a p-dimensional error vector which is uncorrelated with the latent factor. In what follows, we assume that dg P N is a ﬁxed number. Further, Fg “ pfg1 , . . . , fgp qJ denotes a loading matrix where for each j P t1, . . . , pu,
fgj “ pfgj1 , . . . , fgjdg qJ P Rdg is a non-random vector, and Ψg “ diagpψg1 , . . . , ψgp q is a
non-random p ˆ p diagonal matrix whose elements are ψg1 ą 0, . . . , ψgp ą 0. For thelatent vector zgi and error vector gi , we further assume that zgi are iid with Epzgi q “ 0,Epz2gi q “ 1 and Epz4gi q “ κzg ă 8, and gij are iid with Epgij q “ 0, Ep2gij q “ 1 andEp4gij q “ κg ă 8 for g P t1, . . . , ku, i P t1, . . . , ng u, j P t1, . . . , pu and  P t1, . . . , dg u.
Structural assumptions of the model (1) imply that
Epxgi q “ μg , covpxgi q “ Fg FJ
g ` Ψg :“ Σg ,

(2)

pˆp
where Σg P Rpˆp
ą0 and Rą0 denotes the space of real, symmetric, positive deﬁnite, p ˆ pmatrices.
By using the data generated by (1), we design a high-dimensinal test procedure for ageneral linear hypothesis testing (GLHT) problem:rr
H : GM“ O, A : GM
‰ O,

(3)

r is a q ˆ k known coeﬃcient matrix
where M “ pμ1 , . . . , μk qJ is a k ˆ p matrix and G
r to be any pk ´ 1q ˆ k contrast matrix, i.e., anywith full row rank q ă k. By setting G
pk ´ 1q ˆ k matrix with linearly independent rows and zero row sums, the GLHT problem(3) reduces to the one-way MANOVA problem:
H : μ1 “ ¨ ¨ ¨ “ μk , A :‰ H.
Also, various post hoc and contrast tests can be written in the form of (3).

(4)

From Zhang et al. (2017) and Zhang et al. (2022), we re-write (3) into the following
equivalent form:
H : Cμ “ 0, A : Cμ ‰ 0,

(5)

r with D “ diagp1{n1 , ¨ ¨ ¨ , 1{nk q
r G
r J q´1{2 G
where C “ G b Ip (qp ˆ kp matrix), G “ pGD
J J
and μ “ pμJ
1 , . . . , μk q .
řn gJ J
p “ pxJ
Let H “ GJ G and μ
1 , . . . , xk q where xg “ p1{ng qi“1 xgi for g P t1, . . . , ku.
Then, tor testing (5), we deﬁned the test statistic as
#
+
k
ÿ
1
{
}Cp
μ}2 ´
Tnh “
agg trpΨ
gq
p
g“1
#
+
k
ÿ
1
{
p J pH b Ip qpμμ´
agg trpΨ“gq
pg“1
where, for g P t1, . . . , ku, agg is the diagonal element of the matrix A “ D1{2 HD1{2řdpg{and trpΨ
λ pS q. Here, λ pSg q is the th largest eigenvalue of matrix
g q “ trpSg q ´řng “1  g
pxgi ´ xg qpxgi ´ xg qJ and dpg is a consistent estimator of dg based
Sg “ t1{png ´ 1qui“1on the ER method proposed by Ahn and Horenstein (2013). Besides, we derived thelimiting null distribution of Tnh under some assumptions and constructed test procedure
for testing p5q. Also, we compared, through simulations, the performance of the proposed
test and existing procedures suitable for one-way MANOVA problem in high-dimensional
data in terms of size control and power.
References[1] Ahn, S. C., Horenstein, A R., 2013. Eigenvalue ratio test for the number of factors. Econometrica, 81, 1203–1227.
[2] Zhang, J.-T., Guo, J., Zhou, B., 2017. Linear hypothesis testing in high-dimensional oneway MANOVA. J. Multivar. Anal., 155, 200–216.
[3] Zhang, J.-T., Zhou, B., Guo, J., 2022. Linear hypothesis testing in high-dimensional heteroscedastic one-way MANOVA: A normal reference L2 -norm based test. J. Multivar.
Anal., 187, 104816.