Osaka Central Advanced Mathematical Institute (OCAMI)
Osaka Metropolitan University
MEXT Joint Usage/Research Center on Mathematics and Theoretical Physics

OCAMI Reports Vol. 8 (2022)
doi: 10.24544/ocu.20221208-007

Mathematical optimization and statisticaltheories using geometric methods
Organized byHideto Nakashima
Yoshihiko Konno
Hideyuki Ishi
Kenji Fukumizu
October 20–21, 2022
Abstract
This workshop was held on October 20–21, 2022 in order to connect researchers
in several fields, in particular Statistics, Machine Learning and Mathematics, andto share problems and researches in these fields interdisciplinary.

2020 Mathematics Subject Classification.
20G05, 22F30, 43A85, 60E05, 62E10,
62H12, 62J05, 62J07, 62R01
Key words and Phrases.
Algebraic statistics, LASSO, SLOPE, exponential families,
machine learning, geometric analysis© 2022 OCAMI.
OCAMI. Mathematical optimization and statistical theories using geometric methods. OCAMI Reports.
Vol. 8, Osaka Central Advanced Mathematical Institute, Osaka Metropolitan University. 2022, 165 pp.
doi: 10.24544/ocu.20221208-007

Preface
This is a proceedings of the international workshop “Mathematical optimization andstatistical theories using geometric methods” held from October 20th to October 21st in
2022. This workshop aimed to connect researchers in several fields, in particular Statistics,
Machine Learning and Mathematics, and to share problems and researches in these fields
interdisciplinary.
This workshop was supported by Osaka Metropolitan University, Advanced Mathematical Institute MEXT Joint Usage/Research Center on Mathematics and Theoretical
Physics, and also supported by Japan Science and Technology Agency, CREST: “Innovation of Deep Structured Models with Representation of Mathematical Intelligence” in“Creating information utilization platform by integrating mathematical and information
sciences, and development to society.”
This workshop was held in a hybrid format. Domestic speakers are gathered in Academic Extension Center (Osaka Metropolitan University), Foreign speakers participatedby Zoom. We had 10 talks, 6 of which were from Japan and the others were from abroad,and 26 people had been registered in this workshop.

Organizers
Hideto Nakashima
Research Center for Statistical Machine Learning, The Institute of Statistical Mathematics, 10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japan
Email address: hideto@ism.ac.jp
Yoshihiko Konno
Department of Mathematics, Osaka Metropolitan University, 1-1, Gakunen-cho, Naka-ku,
Sakai-shi, 599-8531
Email address: konno@omu.ac.jp
Hideyuki Ishi
Department of Mathematics, Osaka Metropolitan University, 3-3-138, Sugimoto, Sumiyoshiku, Osaka, 558-8585, Japan
Email address: hideyuki-ishi@omu.ac.jp
Kenji Fukumizu
Research Center for Statistical Machine Learning, The Institute of Statistical Mathematics, 10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japan
Email address: fukumizu@ism.ac.jp

Contents
Shoji Toyota
Invariance Learning based on Label Hierarchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Sho Sonoda
Ridgelet Transforms for Neural Networks on Manifolds and Hilbert Spaces . . . . . . 14
Tomonari Sei and Ushio Tanaka
Stein-type distributions on Riemannian manifolds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
Tomasz Skalski
On LASSO and SLOPE estimators and their pattern recovery . . . . . . . . . . . . . . . . . . 49
Carlos Améndola
Likelihood geometry of correlation models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
Piotr Zwiernik
Mixed convex exponential families and locally associated graphical models . . . . . . . 75
Koichi Tojo
Classification problem of invariant q-exponential families on homogeneous spaces85
Yoshihiko Konno
Adaptive shrinkage of singular values for a low-rank matrix mean when a covariancematrix is unknown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
Satoshi Kuriki
Expected Euler characteristic heuristic for smooth Gaussian random fields with inhomogeneous marginals1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
Piotr Graczyk
Pattern recovery by SLOPE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .123
Program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

1

Slides are not included in this report.

Mathematical Optimization and Statistical Theories using Geometric Methods

Invariance Learning based on Label Hierarchy
Shoji Toyota
The Graduate University for Advanced Studies (SOKENDAI)
Training data used in machine learning may contain features that are spuriously correlated to the labels of data. Deep Neural Networks (DNNs) often learn such biased
correlations embedded in training data and hence may fail to predict desired labels of test
data generated by a diﬀerent distribution from one to provide training data. To solve theproblem, Invariance Learning (IL) is a rapidly developed approach to overcome the issueof biased correlation, which is caused by some bias in the distribution of a training dataset(e.g., [1]). IL estimates a predictor invariant to the change of distributions, aiming at
keeping good performance in unseen distributions as well as in the training distributions.
While the IL approach has attracted much attention, requiring training data frommultiple distributions may hinder wide applications in practice; preparing training data
in many distributions often involves expensive data annotation.
To mitigate the problem of annotation cost, we propose a novel IL framework for thesituation where the training data of target classiﬁcation is given in only one distribution,
while the task of higher label hierarchy, which needs lower annotation cost, has datafrom multiple distributions. The new IL framework signiﬁcantly reduces the annotation
cost in comparison with previous IL methods; we need exhausting annotation of originalclasses only for one distribution and just causaer labels for other distributions. Numericalsimulations and theoretical analysis verify the eﬀectiveness of our framework.

References[1] M. Arjovsky, L. Bottou, I. Gulrajani, and D. Lopez-Paz. Invariant Risk Minimization.
arXiv:1907.02893, 2019.

1

2

OCAMI Reports Vol. 8 (2022)

Invariance Learning
based on Label Hierarchy
Shoji Toyota

The Graduate University for Advanced Studies(Joint work with Prof. Kenji Fukumizu )OCAMI workshop, 20 ~ 21, October, 2022
※ The presentation is based on https://arxiv.org/abs/2203.15549. To appear in Neurips 2022.

AgendaBackground
Mathemathical Formulation
MethodTheoryExperimentMathematical Optimization and Statistical Theories using Geometric Methods

AgendaBackground
Mathemathical Formulation
MethodTheoryExperiment

Recent Problem in Machine Learning :Estimators inherit spurious correlation in training data
[S. Beery et al. 2018]

Train

Test

＋

Teacher label

There are no cows !

3

4

OCAMI Reports Vol. 8 (2022)

Statistical Invariance [Arjovsky et al. 2019, Peters et al. 2016]
Notations: Image,: Label,: Test dist.

: Training dist.,

E

Φ(X)

Y

Random variable designating dist.s

Featured Images

Labels(）

Estimating a feature nap Φ by traing data from multiple training dist.s e1, ..., en

Annotation cost problem in Invariance Estimation
Teacher labels are not often attached in images.

・・・

・・・

User
Teacher label

Teacher label

Annotation Vendor
Teacher label

Bird mBird 1

Cost is high especially when
the number of class is large.

Turtle 1

Turtle n・・・※ Images are Cited from [Wah et al., 2011].

Mathematical Optimization and Statistical Theories using Geometric Methods

Invariance Estimation Based on Label Hierarchy
Proposed Framework: Invariance estimation with the following two data
1aTeacher labela2

Label in higher level

Bird 1a・・・

Bird mLabel in higher level

Turtle n

Turtle 1

・・・

・・・

Bird

TurtleInvariance Learning Based on Label Hierarchy

・・・

・・・

User

a

a

Teacher labelaLabel in higher level

Annotation Vendor

Label in higher level

lower cost !
Labels in higher level※ Note that Z = g(Y) holds for some surjective function g.
※ The relation gives effective estimation method.

asOriginal label

56

OCAMI Reports Vol. 8 (2022)

AgendaBackground
Mathemathical Formulation
MethodTheoryExperiment

Mathemathical Formulation: Index designating a dist.
Available
samples: Image and label on e.

1

2(=g(Yie) ): labels in higher hierarchy

Goal: out-of-distribution (o.o.d.) risk Minimization

Out-of-distribution (o.o.d.) risk
: Risk on
Assumption:

Mathematical Optimization and Statistical Theories using Geometric Methods

AgendaBackground
Mathemathical Formulation
MethodTheoryExperiment

Estimation
Esrimaiton object:1.

Feature map Φ which satisfies E → Φ(X) → Y
We can not estimate it by data on a single domain...

2. classifier w predicting a label Y from Φ(X)7

8

OCAMI Reports Vol. 8 (2022)

Estimation
Esrimaiton object:1.

Feature map Φ which satisfies E → Φ(X) → Y

ZLabels in higher level

2. classifier w predicting a label Y from Φ(X)Method: We estimate Φ and w stemiously,
by minimizing the following objective function.

(Dependence measure of E → Φ(X) → Z)
estimating w : evaluated by original label data

estimating Φ: evaluated by higher label data

※ second term:[M. Arjovsky et al. 2019].

Difficulty of Hyperparameter selection[Galrajani et al. 2021]

If we select λ by a naive CV method using training data,
famous methods result in random guess classifiers....
d ^D[-k] , D[k] : Training and Validation,fλ ← fλ・Cross-Validation (CV) for minimizing an o.o.d. risk maxe R(f)

Goal:
fλ

D[-k]

^
fλ

argminλ maxe Re (fλ)
maxeRe(f)D[k]^maxeRe(f)❌

Difficulity: o.o.d. risk estimation from validation dataMathematical Optimization and Statistical Theories using Geometric Methods

Proposed CV methods
Goal:

How can we estimate a risk on? ( ※ De = { (x, z) } )Method 1：Using a risk w.r.t. higher label data Z alternatively.
De[k]= { (x, z) }Proposed
CVMethods
methods
ProposedCV
How can we estimate risk on
Method 1：Using a risk w.r.t. higher label data Z alternatively.
De[k]Method 2：Risk correction(output: probability, loss: cross-entropy)

Thm. (Decomposition formula of risk)

De1[k] ={(x, y)}

?

9

10

OCAMI Reports Vol. 8 (2022)

AgendaBackground
Mathemathical Formulation
MethodTheoryExperimentThoretical analysis of CV methodsd ^fλ ← fλ ( ※ There are some open problems. ): Approximations of an o.o.d. risk by Method I and II (ignoring estimation).

()

Hyperparameterselected by method I

Optimal hyperparameter

Hyperparameterselected by method II

Optimal hyperparameter

The inclusions represent
the success of CVs(with ignoring estimations ).

Mathematical Optimization and Statistical Theories using Geometric Methods

Thoretical analysis of CV methods: all m'ble funct.
: [Rojas-Carulla et al. 2018 ]Correctness of Method 1 (Simplified)Correctness of Method 1 (Simplified): Method II is more applicable !

AgendaBackground
Mathemathical Formulation
MethodTheoryExperiment

11

12

OCAMI Reports Vol. 8 (2022)

Experiment: Image Recognition with 17 class labels
Modified dataset of BREEDS [S. Santurkar et al. 2021]

o.o.d. benchmark constructed by ImageNet [J. Deng et al. 2009]

Dataset e1

Dataset e2

Birds

1
Animals

17 class

Turtles・・・2
2 class (Animals or Non-animals)

Non-animals

Structures

Our method is validated by the worst acc. among e1 and e2

Result (5 runs)
Proposed Methods v.s. Competitors
0.71
0.6

0.56

0.5
0.4
Min{Test Acc on e1,
Test Acc on e2 }

0.36

0.56

|

0.54

0.52

|

0.36

|

0.3
0.2

|

|

0.16

0.1

|

0.06

|

0
ERM

Transfer I Transfer II

Competitors

DA

Ours + CVI Ours + CVII CV comp. I CV comp. I

OracleMathematical Optimization and Statistical Theories using Geometric Methods

References
1. S. Beery et al. Recognition in terra in cognita. In CCV, 2018.
2. M. Arjovsky et al., Invariant Risk Minimization, arXiv:1907.02893, 2019.
3. J. Peters et al. Causal inference using invariant prediction: identificationand confidence intervals, JRSS-B, 2016.
4. M. Rojas-Carulla et al. Invariant models for causal transfer learning. JMLR, 2018
5. S. Santurkar et al. Breeds: Benchmarks for subpopulation shift. In ICLR, 2021.
6. J. Deng et al. Imagenet: A large-scale hierarchical image database. In CVPR, 2009.

13

14

OCAMI Reports Vol. 8 (2022)

Ridgelet Transforms for Neural Networks on
Manifolds and Hilbert Spaces
Sho Sonoda
RIKEN AIP, Tokyo 103–0027 Japan
sho.sonoda@riken.jp

Abstract
To investigate how neural network parameters are organized and arranged, it iseasier to study the distribution of parameters than to study the parameters in eachneuron. The ridgelet transform is a pseudo-inverse operator (or an analysis operator)that maps a given function f to the parameter distribution γ so that a network
Z
S[γ](x) :=γ(a, b)σ(a · x − b)dadb, x ∈ Rm
Rm ×R

represents f , i.e., S[γ] = f . For depth-2 fully-connected networks on Euclidean space,
the ridgelet transform has been discovered up to the closed-form expression, thus we
could describe how the parameters are organized. However, for a variety of modern
neural network architectures, the closed-form expression has not been known . Recently,
our research group has developed a systematic scheme to derive ridgelet transforms
for fully-connected layers on manifolds (noncompact symmetric spaces G/K) (Sonoda
et al., 2022b) and for group convolution layers on abstract Hilbert spaces H (Sonoda
et al., 2022a). In this talk, the speaker will explain a natural way to derive those ridgelet
transforms.

References
S. Sonoda, I. Ishikawa, and M. Ikeda. Universality of Group Convolutional Neural Networks Basedon Ridgelet Analysis on Groups. In Advances in Neural Information Processing Systems 35, 2022a.
S. Sonoda, I. Ishikawa, and M. Ikeda. Fully-Connected Network on Noncompact Symmetric Space andRidgelet Transform based on Helgason-Fourier Analysis. In Proceedings of the 39th International
Conference on Machine Learning, volume 162, 2022b.

1

Mathematical Optimization and Statistical Theories using Geometric Methods

15

The Ridgelet Transforms of Neural Networks on
Manifolds and Hilbert Spaces
Sho Sonoda
Research Scientist
RIKEN Center for Advanced Intelligence Project (AIP), Tokyo, Japan

Mathematical Optimization and Statistical Theories Using Geometric Methods
Osaka Metropolitan University
October 20-21, 2022

1 / 20Q. What is a typical solution obtained by deep learning?
X

Z1

Z2

Z3

Z4

Z5

Y

“camel, 98%”
“rug, 88%”

RGB valuesx ∈ Rminformation processing

label prob.
y ∈ RkWant to identify what solution is typically acquired via deep learning
Want to know why (and when) deep learning performs better (than shallow networks)

2 / 20

16

OCAMI Reports Vol. 8 (2022)

Reparametrization
Finite-width (Discrete, or “Ordinary”) NN

SNN(x; θd ) =

d
X
i=1

ci σ(ai · x − bi )nonlinear parameters: θd = {(ai , bi , ci )}di=1 ∈ R(m+2)dcontinuum limit
Infinite-width (Continuous, or Integral Representation of) NN

S[γ](x) =

Z

Rm ×R

discretizationPd
γd = i=1 ci δ(ai ,bi )γ(a, b)σ(a · x − b)dadb

linear parameter: γ ∈ Map(Rm × R → C)
3 / 20Definition (Ridgelet Transform)
For any function f : Rm → C and ρ : R → C, put
ZR[f ; ρ](a, b) =f (x)ρ(a · x − b)dx,
Rm(a, b) ∈ Rm × R.

Theorem (Reconstruction Formula)
For any σ ∈ S 0 (R), ρ ∈ S(R) and f ∈ L2 (Rm ), we haveZ
S[R[f ; ρ]](x) = R[f ; ρ](a, b)σ(a · x − b)dadb = ((σ, ρ))f (x),R
where ((σ, ρ)) = (2π)m−1 R σ ] (ω)ρ] (ω)|ω|−m dω and ] denotes the Fourier transform

Meaning 1: Continuous NN is a universal approximator
Meaning 2: R and S play the same role as Fourier F and inverse Fourier F −1 transforms:
Z−1−m
F [F [f ]](x) = (2π)
F [f ](ξ)eix·ξ dξ = f (x)Rm

Independently “discovered” by Murata (1996), Candès (1998), and Rubin (1998)

4 / 20Mathematical Optimization and Statistical Theories using Geometric Methods

17

Numerical example of ridgelet transform R[f ; ρ](a, b)

−1.0

−0.5

y
0.0

0.5

1.0

f (x) = sin(2πx)
R 1[−1,1] (x)P
R[f ; ρ](a, b) = R f (x)ρ(ax − b)dx ≈ i sin(2πxi )ρ(axi − b)∆x
σ(b) = tanh(b)(2)
ρ(b) = H[ρ0 ](b) with ρ0 (b) := exp(−b2 /2), Hilbert transform H

−2

−1

0
x

1

2

data f (x)image R[f ; ρ](a, b)

5 / 20

Visualization results of reconstruction formula S[R[f ; ρ]] = ((σ, ρ))f(1)

ρ1 = Hρ0((σ, ρ1 )) = 0(2)

ρ2 = Hρ0((σ, ρ2 )) 6= 0

(3)

ρ3 = Hρ0((σ, ρ3 )) = 0(4)

ρ4 = Hρ0((σ, ρ4 )) 6= 0

6 / 20

18

OCAMI Reports Vol. 8 (2022)

How the parameter distribution looks like?

−1.0

−0.5

y
0.0

0.5

1.0

Pd
We will train many (n = 1, 000) neural networks SNN(x; θd ) = j=1 cj σ(aj · x − bj ) with d = 10
hidden units, and see the distribution of trained parameters (aj , bj , cj ).
Data generating function: f (x) = sin(2πx) 1[−1,1] (x)σ(z) = tanh(z)SGD w. weight decay

−2

−1

0
x

1

2

data f (x)7 / 20A scatter plot of d × n = 10 hidden parameters (aj , bj , cj ) obtained from n = 1 neural networkPd
j=1 cj σ(aj · x − bj ) with d = 10 hidden units.

b n (θd )arg min L

7 / 20Mathematical Optimization and Statistical Theories using Geometric Methods

19

A scatter plot of d × n = 20 hidden parameters (aj , bj , cj ) obtained from n = 2 neural networks withd = 10 hidden units.

b n (θd )arg min L

7 / 20A scatter plot of d × n = 50 hidden parameters (aj , bj , cj ) obtained from n = 5 neural networks withd = 10 hidden units.

b n (θd )arg min L

7 / 20

20

OCAMI Reports Vol. 8 (2022)

A scatter plot of d × n = 100 hidden parameters (aj , bj , cj ) obtained from n = 10 neural networks withd = 10 hidden units.

b n (θd )arg min L

7 / 20A scatter plot of d × n = 500 hidden parameters (aj , bj , cj ) obtained from n = 50 neural networks withd = 10 hidden units.

b n (θd )arg min L

7 / 20Mathematical Optimization and Statistical Theories using Geometric Methods

21

A scatter plot of d × n = 1, 000 hidden parameters (aj , bj , cj ) obtained from n = 100 neural networks
with d = 10 hidden units.

×100
b n (θd )arg min L

7 / 20A scatter plot of d × n = 5, 000 hidden parameters (aj , bj , cj ) obtained from n = 500 neural networks
with d = 10 hidden units.

×500
b n (θd )arg min L

7 / 20

22

OCAMI Reports Vol. 8 (2022)

A scatter plot of d × n = 10, 000 hidden parameters (aj , bj , cj ) obtained from n = 1, 000 neuralnetworks with d = 10 hidden units.

×1000
b n (θd )arg min L

7 / 20appears to be the image R[f ; ρ] of data f .
(∞)
(0) b
(formal) θ:= SGD(θ , Ln ) ∼ R[f ; ρ] (including sign!)ddb n (θd )scatter plot arg min L

image R[f ](a, b)
7 / 20Mathematical Optimization and Statistical Theories using Geometric Methods

23

Q. How to Find R?—A. Solve S[γ] = f
Appendix A.3, in Sonoda-Ishikawa-Ikeda, arXiv:2106.04770

Step 1. Turn the network into a Fourier expressionZ Z
S[γ](x) =γ(a, b)σ(a · x−b)db da
Rm
R

Z
Z 
R ]
1
1
γ ] (a, ω)σ ] (ω)eiωa·x dω da, ∵ 2π
γ (a, ω)σ ] (ω)eiωb dω = (γ(a, •) ∗ σ)(b)=R
Rm 2π RZ Z
1=γ ] (ξ/ω, ω)eiξ·x dξ |ω|−m σ ] (ω)dω, by (a, ω) = (ξ/ω, ω)
2π R Rm
where ·] is the Fourier transform in b
Step 2. Assume a separation-of-variables form]
γf,ρ(ξ/ω, ω) := fb(ξ)ρ] (ω)

Then, (1) γf,ρ is a particular solution
Z
 Z

1
]
−m
iξ·x
]
b
S[γf,ρ ] =σ (ω)ρ (ω)|ω| dω
f (ξ)edξ = ((σ, ρ))f2π(2) and γf,ρ (a, b) = R[f ; ρ](a, b).

8 / 20Further Results
Theorem (S-Ishikawa-Ikeda, AISTATS2021)
The empirical regularized least squares parameters in the finite NNs converges to the ridgelet transform:n1Xn,d→∞,β→+0|f (xi ) − S[γd ](xi )|2 + β|c|2 −−−−−−−−−→ S ∗ [f ] = R[f ; σ∗ ]nci δ(a ,b )i=1

arg min
P

γd =

d
i=1

i

i

Ridgelet transform can characterize the parameters obtained by learning (loss minimization)

Theorem (S-Ishikawa-Ikeda, arXiv:2106.04770)
The general solution of S[γ] = f is given by a sum of ridgelet transforms
X
γ = S ∗ [f ] +cij R[ei ; ρj ]ijwhere ei and ρj are ONSs in L2 (Rm ) and L2 (R, ((·, ·))) resp. satisfying ((σ, ρj )) = 0
Ridgelet transform is not only sufficient but also necessary
9 / 20

24

OCAMI Reports Vol. 8 (2022)

Extensions to modern network architectures

Based on the Fourier expression technique, we have developed new ridgelet transforms for
1

Group convolutional NNs on Hilbert space H
in S-Ishikawa-Ikeda (NeurIPS2022) and2

Fully-connected NNs on manifold (noncompact symmetric space) G/K
in S-Ishikawa-Ikeda (ICML2022)

10 / 20

Group Convolutional NNs on Hilbert Space H1

1 S-Ishikawa-Ikeda, NeurIPS2022
11 / 20

Mathematical Optimization and Statistical Theories using Geometric Methods

25

Definition (Group CNN)
Let G be a group, H be a Hilbert space, and T : G → GL(H) be a group representation. Let Hm ⊂ H
be an m-dimensional subspace equipped with the Lebesgue measure λ. PutZ
S[γ](x)(g) :=γ(a, b)σ((a ∗ x)(g) − b)dλ(a)db, x ∈ H, g ∈ G
Hm ×R

where the (G, T )-convolution is given by(a ∗ x)(g) := hTg−1 [x], aiH .

Example (Cyclic CNN for multichannel image)0

CNN(x)(p, q) =

n
X
`=1



c` σ 

n X
mX

k=1 i,j=1

2



k
`
ak`
,
ij xi+p,j+q − b2

x = (xkij ) ∈ Rm ×n , (p, q) ∈ (Z/mZ)2

i.e., G = (Z/mZ)2 , H = Rm ×n , Tp,q (x) := (x••−p,•−q )12 / 20In the following, e ∈ G denotes the identity element.

Definition (Ridgelet Transform)
For any function f : Hm → CG and ρ : R → C, put
ZR[f ; ρ](a, b) :=f (x)(e)ρ(ha, xiH − b)dλ(x).
HmDefinition ((G, T )-Equivariance)
A (nonlinear) map f : H → CG is (G, T )-equivariant whenf (Tg [x])(h) = f (x)(g −1 h),∀x ∈ Hm , g, h ∈ GTheorem (Reconstruction Formula)
Suppose that f is (G, T )-equivariant and f (•)(e) ∈ L2 (Hm ), then S[R[f ; ρ]] = ((σ, ρ))f .
Meaning: Universality of continuous GCNN
Corollary: cc-universality of finite GCNNs
13 / 20

26

OCAMI Reports Vol. 8 (2022)

Sketch Proof
Step 1. Turn to Fourier expression:Z
S[γ](x)(g) =γ(a, b)σ(hTg−1 [x], aiH − b)dadb
Hm ×R
Z
1
=
γ ] (a, ω)σ ] (ω)eiωhTg−1 [x],aiH dadω
2π Hm ×RZ
1=γ ] (ξ/ω, ω)σ ] (ω)eihTg−1 [x],ξiH |ω|−m dξdω.
2π Hm ×R
Step 2. Put separation-of-variables form:]
γf,ρ
(ξ/ω, ω) := fb(ξ)(e)ρ] (ω).

By the construction it is a particular solution:
Z
Z
1
S[γf,ρ ](x)(g) =fb(ξ)(e)eihTg−1 [x],ξiH dλ(ξ) σ ] (ω)ρ] (ω)|ω|−m dω
2π HmR= ((σ, ρ))f (x)(g).

and γf,ρ = R[f ; ρ].
14 / 20

Fully-Connected NNs on Noncompact Symmetric Space2

2 S-Ishikawa-Ikeda, ICML2022
15 / 20

Mathematical Optimization and Statistical Theories using Geometric Methods

27

Definition (Fully-Connected NNs on Noncompact Symmetric Space G/K)
Let G be a connected semisimple real Lie group, let G = KAN be the Iwasawa decomposition, and let
X := G/K be the noncompact symmetric space. PutZ
S[γ](x) :=γ(a, u, b)σ(ahx, ui − b)e%hx,ui dadudb, x ∈ X = G/Ka∗ ×∂X×R

where a∗ is the dual of Lie algebra of A, ∂X is the boundary, and hx, ui is an X-counter of the
Euclidean inner product x · u for (x, u) ∈ Rm × Sm−1 .

Example (Continuous Horospherical Hyperbolic NN)mm
On the Poincaré ballPmmodel B := {x ∈ R | |x| < 1} equipped with the Riemannian metric−2
g = 4(1 − |x|)i=1 dxi ⊗ dxi ,Z
S[γ](x) :=γ(a, u, b)σ(ahx, ui − b)e%hx,ui dadudb, x ∈ Bm
R×∂Bm ×R

% = (m − 1)/2, hx, ui = log



1 − |x|2E|x − u|2E,(x, u) ∈ Bm × ∂Bm
16 / 20

Definition (Ridgelet Transform)
For any function f : X → C and an auxiliary function ρ : R → C, put
ZR[f ; ρ](a, u, b) :=c[f ](x)ρ(ahx, ui − b)e%hx,ui dx
Xwhere c[f ] is a Helgason-Fourier multiplier.

Theorem (Reconstruction Formula)
For any σ ∈ S 0 (R), ρ ∈ S(R), and f ∈ L2 (X), we haveZ
S[R[f ; ρ]] =R[f ; ρ](a, u, b)σ(ahx, ui − b)e%hx,ui dadudb = ((σ, ρ))f.
a∗ ×∂X×R

where ((σ, ρ)) is a certain scalar product.
Meaning: Universality of continuous Fully-Connected NN on X
Corollary: cc-universality of finite Fully-Connected NNs on X
17 / 20

28

OCAMI Reports Vol. 8 (2022)

Fourier Analysis on X = G/KHelgason, GGA (1984, Introduction); GASS (2008, Chapter III)Definition (Helgason-Fourier Transform)
For any function f : X → C,
fb(λ, u) :=Z

f (x)e(−iλ+%)hx,ui dx,

X(λ, u) ∈ a∗ × ∂X

with a certain constant vector % ∈ a∗ .

Theorem (Inversion Formula)
For any f ∈ L2 (X) (or f ∈ Cc∞ (X)),
Z
f (x) = |W |−1a∗ ×∂Xfb(λ, u)e(iλ+%)hx,ui |c(λ)|−2 dλdu,x∈X

where c is the Harish-Chandra c-function, and |W | is a constant.

This is a “Fourier transform” because e(−iλ+%)hx,ui is the eigenfunction e(−iλ+%)hx,ui of the
Laplace-Beltrami operator ∆X on X

18 / 20

Sketch Proof
Given a function f : G/K → C, consider solving an integral equation S[γ] = f of unknown γ.
Step 1: Change the frame of S[γ] from neurons to a Fourier expression:Z
S[γ](x) :=γ(a, u, b)σ(ahx, ui − b)e%hx,ui dadudba∗ ×∂X×RZ Z
1]2 (iλ+%)hx,ui dλdu=γ (λ/ω, u, ω)|c(λ)| e|ω|−r σ ] (ω)dω,
2π R a∗ ×∂X|c(λ)|2
where ] denotes the Euclidean-Fourier transform in b.
Step 2: Since inside [· · · ] is the inverse Helgason-Fourier transform, put a separation-of-variables
form:
]
γf,ρ
(λ/ω, u, ω) = fb(λ, u)ρ] (ω)|c(λ)|−2 .

Then, by the construction, it is a particular solution:

S[γf,ρ ] = ((σ, ρ))f,|W | R

where ((σ, ρ)) := 2π R σ ] (ω)ρ] (ω)|ω|−m dω.
In the end, we can verify that γf,ρ is the ridgelet transform R[f ; ρ].
19 / 20

Mathematical Optimization and Statistical Theories using Geometric Methods

29

Conclusion
Ultimate goal:
I

Characterize deep solutions

We have seen:I

Shallow solutions are characterized by ridgelet transform

Take home message:
I

If there is a Fourier transform, then so is the ridgelet transform

We will see:I

A ridgelet transform for depth

20 / 20

30

OCAMI Reports Vol. 8 (2022)

Stein-type distributions on Riemannian manifolds
Tomonari Sei (The University of Tokyo)∗1
Ushio Tanaka (Osaka Metropolitan University)∗2

1. Stein-type distributions on the Euclidean space
Let P 2 be the set of probability distributions µ on Rd with mean zero and finite second
moments such that each marginal distribution µi (i = 1, . . . , d) is absolutely continuouswith respect to the Lebesgue measure dxi on R. We say that a probability distributionµ ∈ P 2 is Stein-type if it satisfies
( d
)
∫
∫
∑
f (xi )
xj dµ = f ′ (xi )dµ, i = 1, . . . , d,j=1

for any absolutely continuous function f : R → R with bounded derivative f ′ .
Let Tcw be the set of coordinate-wise transformations T (x) = (T1 (x1 ), . . . , Td (xd ))such that each Ti is non-decreasing. In [2], it is shown that for any given µ0 ∈ P 2 ,there exists T ∈ Tcw such that T♯ µ0 is Stein-type. The transformation is characterized
by a minimizer of a functional
)2
∫ (∑
d ∫
d
∑
dµi
1
F (µ) =logdµi +xi dµ,
dxi
2 i=1
i=1
over a fiber {T♯ µ0 | T ∈ Tcw }. The fiber is totally geodesic in the L2 -Wasserstein spaceand F is convex with respect to displacement interpolation. The optimal map T isapplied to the problem of determining a general index in [2].

2. Generalization to manifolds
We generalize the Stein-type distributions to those∏don Riemannian manifolds. Thedspace R is replaced with a product space M = i=1 Mi , where each Mi is a Riemannian manifold. The space P 2 of distributions is defined as well. Let Tcw be theset of coordinate-wise transformations T (x) = (T1 (x1 ), . . . , Td (xd )) such that each
Ti : Mi → Mi is monotone. Here, Ti is said to be monotone if it is written as
Ti (xi ) = expxi ∇ϕi (xi ) with a cost convex function ϕi : Mi → R (see [1]). The Steintype distribution is defined by a minimizer of a functional
∫
d ∫
∑
dµi
F (µ) =log
dµi + V (x)dµ,dxii=1
over a fiber {T♯ µ0 | T ∈ Tcw }, where V : M → R is a given function.

References[1] McCann, R. J. (2001). Polar factorization of maps on Riemannian manifolds, Geometricand Functional Analysis, 11, 589–608.
[2] Sei, T. (2022). Coordinate-wise transformation of probability distributions to achieve a
Stein-type identity, Information Geometry, 5, 325–354.
∗1∗2e-mail: sei@mist.i.u-tokyo.ac.jp
e-mail: utanaka@omu.ac.jp

Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

31
Summary

Stein-type distributions on Riemannian manifolds
Tomonari SEIThe University of Tokyo

Ushio Tanaka
Osaka Metropolitan University

Oct 20 (Thu), OCAMI workshop“Mathematical optimization and statistical theories
using geometric methods”

1 / 35Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

The Stein identity
We begin with the following fact.
Proposition (Stein identity)
A random variable X follows N(0, 1) if and only if
E [Xf (X )] = E [f ′ (X )]for any diﬀerentiable function f with bounded f ′ .
2

Proof: (⇒) For the density function ϕ(x) = (2π)−1/2 e −x /2 ,∫∫∫′xf (x)ϕ(x)dx = f (x){−ϕ(x)} dx = f ′ (x)ϕ(x)dx.

(⇐) If E [Xf (X )] = E [f ′ (X )], it is shown that X has density p(x).
Then the identity is equivalent top ′ (x) + xp(x) = 0.
The unique solution is p(x) = ϕ(x).

2 / 35

32

OCAMI Reports Vol. 8 (2022)

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Application of the Stein identity
Why is the Stein identity important?
Stein’s unbiased risk estimator (statistics)
Central limit theorem (probability theory)
Stein discrepancy (machine learning)
Application: Stein’s unbiased risk estimator
Let X ∼ Nd (θ, Id ), where θ ∈ Rd is unknown parameter.
Consider an estimator X + f (X ) of θ. The risk is
E [∥X + f (x) − θ∥2 ]= E [∥X − θ∥2 ] + 2E [f (X )⊤ (X − θ)] + E [∥f (X )∥2 ]= d + 2E [∇⊤ f (X )] + E [∥f (X )∥2 ] (Stein identity)= E [d + 2∇⊤ f (X ) + ∥f (X )∥2 ]|{z}risk estimator
Introduction

OGI

Stein-type distribution

Known results

3 / 35Generalization to manifolds

Summary

Another application of Stein identity
Application: Poincaré inequality (Chernoﬀ 1981, Chen 1982)
If X ∼ N(0, 1), thenV [g (X )] ≤ E [g ′ (X )2 ],with equality if and only if g (x) = ax + b.
Indeed,V [g (X )] ≤ E [(g (X ) − g (0))2 ]∫X= E [( 0 g ′ (x)dx)2 ]∫X≤ E [X 0 g ′ (x)2 dx] (Cauchy–Schwarz*)= E [g ′ (X )2 ] (Stein identity).

(* valid even for X < 0.)
4 / 35Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

33
Summary

Outline of this talkIn this talk, we generalize the Stein identity in the following
manner.
1

Define Stein-type distributions on Rd by an identity
E [(X1 + · · · + Xd )f (Xi )] = E [f ′ (Xi )].

2

Define Stein-type distributions on the direct product ofRiemaniann manifolds (on-going work).

We first see the background of the problem in a couple of slides.

5 / 35Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Background: Objective general index (OGI)
S. (2016) pointed out that the Stein identity is related to a scaling
problem, which is a motivation of this work.
First, consider d random variables X1 , . . . , Xd .
For example, Xi is academic score of students on i-th subject.
Proposition (S. 2016)
There exist unique w1 , . . . , wd > 0 such that
Cov(Y , wi Xi ) = 1(i = 1, . . . , d),where Y = w1 X1 + · · · + wd Xd , under a mild condition.
The proof is based on matrix scaling (Marshall–Olkin 1968).
We call Y the objective general index (OGI).
The Stein identity appears in a functional version of this fact.
6 / 35

34

OCAMI Reports Vol. 8 (2022)

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Illustration
A numerical example
Suppose that the covariance matrix of X1 , X2 , X3 is
1−0.5 −0.5
10 .
(Cov(Xi , Xj ))3i,j=1 = −0.5−0.5
01
In this case,

Cov(X1 + X2 + X3 , X1 ) = 1 − 0.5 − 0.5 = 0.
But, a weight (w1 , w2 , w3 ) = (2.135779, 1.667566, 1.667566) gives
Cov(w1 X1 + w2 X2 + w3 X3 , wi Xi ) = 1,
|
{z
}
OGI

Introduction

OGI

Stein-type distribution

Known results

i = 1, 2, 3.

7 / 35Generalization to manifolds

Summary

Functional OGI
Next, consider a random variable X with density p(x).
Define an infinite number of variables by Heaviside function:hξ (X ) = I{X ≥ξ} − E [I{X ≥ξ} ],ξ ∈ R.

What is OGI of {hξ (X )}ξ∈R ?
Proposition (S. 2016)
There exists a unique positive function w (ξ) such that
Cov(Y , w (ξ)hξ (X )) = 1
where Y =∫(ξ ∈ R),(∗)

R w (ξ)hξ (X )p(ξ)dξ. In fact, Y ∼ N(0, 1).

We call Y the functional OGI of X .
The identity (∗) is considered as a version of the Stein identity.
Let us check it.
8 / 35Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

35
Summary

Functional OGI and Stein identity
It is shown that the condition of the functional OGI
Cov(Y , w (ξ)hξ (X )) = 1
is equivalent to the Stein identity
E [Yfξ (Y )] = E [fξ′ (Y )]∫for fξ (y ) = hξ (T −1 (y )) and T (x) = R w (ξ)hξ (x)p(ξ)dξ.

In other words, the functional OGI is characterized by anincreasing function T that attains the Stein identity.

The Stein-type distribution we now discuss is a generalizationof N(0, 1) based on this fact.

9 / 35Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Variational characterization
Before proceeding, we recall a variational characterization ofN(0, 1).
Proposition∫∫

x2
F (p) =
p(x) log p(x)dx +p(x)dxRR 2
2

has a unique minimizer p(x) = ϕ(x) = (2π)−1/2 e −x /2 .
Proof 1: F (p) =∫p(x) log(p(x)/ϕ(x))dx + const

Proof 2: Let p0 be a minimizer of F . Let T (x) = x + εf (x)be an increasing function. Then,( ∫)∫F (T♯ p0 )−F (p0 ) = ε − p0 (x)f ′ (x)dx + f (x)p0 (x)dx +o(ε).
The stationary condition is the Stein identity. So p0 = ϕ.

10 / 35

36

OCAMI Reports Vol. 8 (2022)

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

SummaryFiber
Now let us go on to the Rd case.
Let P 2 be the set of probability distributions µ on Rd suchthat the marginal distribution µi satisfies
∫
∫ 2
x
dµ
=
0,
µi ≪ Leb.
i
R i
R xi dµi < ∞,We call T : Rd → Rd a coordinate-wise transformation if
T (x) = (T1 (x1 ), . . . , Td (xd )),Ti′ (xi ) > 0.

For each µ ∈ P 2 , define the µ-fiber
Fµ = {T♯ µ ∈ P 2 | T is coordinate-wise},
where T♯ denotes the push forward.
11 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

PictureThe space P 2 is decomposed into the set of fibers.
We define a Stein-type distribution in each fiber.
2

Stein-type

copulasfiber

Remark
Fµ is totally geodesic in the Wasserstein space.

Fµ has a unique copula (Sklar’s theorem). A copula refers toa distribution with uniform marginals on [0, 1].
12 / 35

Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

37
Summary

A free-energy functional
Define a functional F : P 2 → R by
F (µ) =

d ∫
∑
i=1

dµi
1
log
dµi +
dxi
2
R

( d
∑

∫

Rd

i=1

Xi

)2

dµ,

We can further consider
F (µ) =

d ∫
∑
i=1

dµi
log
dµi +
dxi
R

∫

Rd

V (x)dµ,with some smooth function V (x) (S. 2017).
This appears in the optimal transport theory (McCann
1997)
∫
dµ
except that the entropy term is replaced with log dx dµ.
13 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Minimization over the fiberF

F is not bounded from below on the whole space P 2 .

But F may be bounded from below on each fiber Fµ .

14 / 35

38

OCAMI Reports Vol. 8 (2022)

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Stein-type distribution
To minimize F over the fiber, consider a perturbation of thetransformation around the identity:Tε = Id + εf ,f (x) = (f1 (x1 ), . . . , fd (xd )),ε ∈ R.

Then we have, as ε → 0, the first variation∑∫
F ((Tε )♯ µ) ≃ F (µ)+ε{−fi ′ (xi )+fi (xi )(x1 +· · ·+xd )}dµ.
i

Definition (Stein-type distribution)
A distribution µ is called a Stein-type distribution if it satisfies∫∫(x1 + · · · + xd )fi (xi )dµ =fi ′ (xi )dµ, ∀i, ∀fi ∈ C1 (R).
Rd

Rd

15 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

ExamplesExample 1 (independent case)
If X1 , . . . , Xd are independent and have zero mean, then theequation
E [(X1 + · · · + Xd )f (Xi )] = E [f ′ (Xi )]forcesE [Xi f (Xi )] = E [f ′ (Xi )].
Thus, only the independent Stein-type distribution is the standard
normal distribution.

16 / 35

Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

39
Summary

ExamplesExample 2 (Gaussian)Let (X1 , . . . , Xd ) ∼ Nd (0, S). Then the distribution is Stein-type if
and only if


d
∑
∑
Sij = Cov Xi ,Xj  = 1j=1jfor i = 1, . . . , d. This is the same as the OGI property.
Example 3 (non-Gaussian)

Let Z ∼ N(0, 1) and U be any distribution with E [U] = 0 andE [U 2 ] < ∞. Then the random vector (X1 , X2 ) with
X1 =

Z +U
√ ,
2

X2 =

Z −U
√
2

is Stein-type.
17 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Functional OGI (revisited)
We briefly discuss an application of our results.
Problem
Let X1 , . . . , Xd be random variables with joint density p(x),which represent students’ scores on d academic subjects.
How to define the overall score?
An answer
Let Y =

∑dj=1 Tj (Xj ), where T (X ) is the Stein-type.

Then the Heaviside function f (xi ) = hξ (xi ) yields
E [Y | Xi > ξ] > E [Y | Xi < ξ],∀ξ ∈ R, ∀i.

Interpretation: students with higher score on each subject i
has higher overall score in mean.
18 / 35

40

OCAMI Reports Vol. 8 (2022)

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Assumption on µ: copositivity
S. (2022) established an existence and uniqueness theorem. We
suppose some conditions.
For each µ ∈ P 2 , denote the product measure of marginal
distributions by∏µ⊥ = di=1 µi .

Definition (Copositivity)

We say that µ is copositive if∫ ∑{Ti (xi )}2 dµ> 0.
β(µ) = inf ∫ ∑ i
T :cw { i Ti (xi )}2 dµ⊥

Trivially, if µ is independent (µ = µ⊥ ), then β(µ) = 1.
Suﬃcient conditions for copositivity are discussed later.
19 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Assumption on µ: regular support
Definition (Regularity)
We say that µ has a regular support if the support of µ is the
direct product of the supports of µi ’s.

regularnon-regular

20 / 35

Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

41
Summary

Existence and uniqueness theorem
Theorem (Existence and uniqueness)
Suppose that µ is copositive and has a regular support. Then there
exists a unique Stein-type distribution in the µ-fiber.

F

21 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Proof sketch.
Uniqueness follows from the displacement convexity
F ([(1 − λ)T0 + λT1 ]♯ µ) > (1 − λ)F ((T0 )♯ µ) + λF ((T1 )♯ µ),where strict inequality follows from the regular supportcondition.
For existence, we use the copositivity to obtain
F (µ) ≥∫dµ⊥ ⊥ β
log
dµ +
dx
2

∫ (∑ )2xi
dµ⊥ .
iThen the problem is essentially reduced to the independent
case µ = µ⊥ .

22 / 35

42

OCAMI Reports Vol. 8 (2022)

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Suﬃcient conditions for copositivity
We establish suﬃcient conditions for copositivity∫ ∑{Ti (xi )}2 dµβ(µ) = inf ∫ ∑ i> 0.
T :cw { i Ti (xi )}2 dµ⊥

The notion of positive dependence plays a significant role.
Definition (e.g. Rüschendorf 2013)1

2

3

p(x) is called MTP2 (multivariate totally positive of order 2)
if p(x ∨ y )p(x ∧ y ) ≥ p(x)p(y ) for all x, y ∈ Rd .
∫∫∫p(x) is said to be associated if ϕψpdx ≥ ϕpdx ψpdx
for all increasing ϕ, ψ : Rd → R.
p(x)is called PSMD
super-modular dependent) if∫∫ (positive⊥
ϕ(x)p(x)dx ≥ ϕ(x)p (x)dx for any super-modular
function ϕ.

23 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Suﬃcient conditions

Theorem (FKG 1971, Christofides 2004, S. 2017)
MTP2 ⇒ associated ⇒ PSMD ⇒ copositive.
MTP2 is relatively easy to confrim.

24 / 35

Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

43
Summary

Open problems

There are some open problems.
Conjectures
1

The marginal support of any Stein-type density is R.

2

Existence implies uniqueness.

3

A Gaussian distribution is copositive if the covariance matrix
is strictly copositive.

For the rest of talk, we generalize the Stein-type distributions on
Rd to the direct space of Riemannian maniolds.

25 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Optimal transport on Riemannian manifolds
We recall the optimal transport theory on Riemannian manifolds
according to McCann (2001).
Let (M, g ) be a Riemannian manifold that is C 3 , compact
and connected without boundaries.
An example in mind is M = S 1 (circle).
Let d(x, y ) be the geodesic distance between x, y ∈ M.
A cost is defined by c(x, y ) = d(x, y )2 /2.

A function ϕ : M → R is called cost-convex if there exists
ϕ∗ : M → R such that
ϕ(x) = sup {−c(x, y ) − ϕ∗ (y )}.
y ∈M

If ϕ is cost-convex, it is Lipschitz and therefore is diﬀerentiable
vol-a.e. (Rademacher’s diﬀerentiability theorem).
26 / 35

44

OCAMI Reports Vol. 8 (2022)

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

McCann’s theorem

For cost-convex ϕ, a map T : M → M defined by
T (x) = expx (∇ϕ(x))
is considered as a generalization of increasing functions on R.
Theorem (McCann 2001)
Let µ ≪ vol and ν be probability measures on M. Then there
exists a unique cost-convex function ϕ (up to additive constants)
such that T (x) = expx (∇ϕ(x)) pushes µ forward∫ to ν. This map
is a unique minimizer of the transportation cost c(x, T (x))dµ.

27 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

SummaryFiber
Let M1 , . . . , Md be C 3 compact Riemannian manifolds.
∏
Consider the product space M = di=1 Mi .

Let P be the set of probability distributions µ on M such thatthe marginal distribution µi satisfies µi ≪ voli .
We call T : M → M a coordinate-wise transformation if
T (x) = (T1 (x1 ), . . . , Td (xd )),Ti = expxi (∇ϕi (xi )),where ϕi is cost-convex.
For each µ ∈ P, define the µ-fiber
Fµ = {T♯ µ ∈ P | T is coordinate-wise},
where T♯ denotes the push forward.
28 / 35

Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

45
Summary

PictureThe space P is decomposed into the set of fibers.
We define a Stein-type distribution in each fiber.
2

Stein-type

copulasfiber

Remark: Sklar’s theorem on manifolds
Fµ has a unique “copula”, which refers to a distribution withuniform marginals on Mi . (cf. circula; Jones et al. (2015))
29 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Stein-type distribution on M

Let V (x) be a smooth function on M =Define a functional F : P → R by
F (µ) =

d ∫
∑
i=1

dµi
dµi +
log
dxi
Mi

∏d

i=1 Mi .

∫V (x)dµ.

M

Definition
A Stein-type distribution on M is defined by a minimizer of F (µ)
over a fiber.
Problem: Existence and uniqueness? → future work..
30 / 35

46

OCAMI Reports Vol. 8 (2022)

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

Stationary condition
To minimize F over the fiber, consider a perturbation of thetransformation around the identity:Tε (x) = expx (εf (x)),

f (x) = (f1 (x1 ), . . . , fd (xd )),ε ∈ R.

Then we have the first variation∑∫
F ((Tε )♯ µ) ≃ F (µ) + ε{−∇i fi (xi ) + fi (xi )∇i V (x)}dµ.
i

Lemma
If µ is Stein-type, then∫∫fi (xi )∇i V (x)dµ =∇i fi (xi )dµ,MM∀i, ∀fi ∈ C1 (Mi ).
31 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

ExampleCircular case
Let M1 = · · · = Md = S 1 .

We use the coordinate xi = (cos θi , sin θi ) ∈ Mi .
Consider a function∑
1 ∑V (x) = {( i cos θi )2 + ( i sin θi )2 }.
2

The derivative is ∂θi V (x) = −A(θ) sin(θi − θ̄), where A(θ)
and θ̄ are defined appropriately.
Then the Stein-type distribution has to satisfy∫∫− M fi (θi )A(θ) sin(θi − θ̄)dµ = M fi ′ (θi )dµ.
Any application? → future work...

32 / 35

Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

47
Summary

Summary and future work
Summary
We defined the Stein-type distributions on Euclidean spaceand established the existence and uniqueness theorem.
We generalized it to distributions on Riemannian manifolds.
Future works
Existence seems OK due to the compactness. Uniqueness may
be non-trivial.
Any analogue of Poincaré inequality?
We are seeking applications.

Thank you for your attention!

33 / 35

Introduction

OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

References I
Bishop, C. M. (2006). Pattern Recognition and Machine Learning, Springer.
Chen, L. H. (1982). An inequality for the multivariate normal distribution, J.
Multivariate Anal., 12, 306–315.
Chen, L. H. Y., Goldstein, L., and Shao, Q. (2011). Normal Approximation by
Stein’s Method, Springer.
Chernoﬀ, H. (1981). A note on an inequality involving the normal distribution,
Ann. Probab., 9 (3), 533–535.
Christofides, T. C. and Vaggelatou, E. (2004). A connection between
supermodular ordering and positive/negative association, J. Multivariate Anal.,
88, 138–151.
Fallat, S., Lauritzen, S., Sadeghi, K., Uhler, C., Wermuth, N., and Zwiernik, P.
(2017). Total positivity in Markov structures, Ann. Statist., 45 (3), 1152–1184.
Fortuin, C. M., Kasteleyn, P. W., and Ginibre, J. (1971). Correlation inequalities
on some partially ordered sets, Comm. Math. Phys., 22, 89–103.
Jones, M. C., Pewsey, A., and Kato, S. (2015). On a class of circulas: copulas
for circular distributions, Ann. Inst. Statist. Math., 67 (5), 843–862.
Marshall, A. W., Olkin, I., (1968). Scaling of matrices to achieve specified row
and column sums. Numer. Math., 12, 83–90.
34 / 35

48
Introduction

OCAMI Reports Vol. 8 (2022)
OGI

Stein-type distribution

Known results

Generalization to manifolds

Summary

References II
McCann, R. J. (1997). A convexity principle for interacting gases, Adv. Math.,
128, 153–179.
McCann, R. J. (2001). Polar factorization of maps on Riemannian manifolds,
Geometric and Functional Analysis, 11, 589–608.
Müller, A. and Stoyan, D. (2002). Comparison Methods for Stochastic Models
and Risks, Wiley.
Nelsen, R. B. (2006). An Introduction to Copulas, 2nd ed., Springer.
Rüschendorf, L. (1981). Characterization of dependence concepts in normal
distributions, Ann. Inst. Statist. Math., 33, 347–359.
Rüschendorf, L. (2013). Mathematical Risk Analysis, Springer.
Sei, T. (2016). An objective general index for multivariate ordered data, J.
Multivariate Anal., 147, 247–264.
Sei, T. (2017). Coordinate-wise transformation and Stein-type densities,
Proceedings of the 3rd Conference on Geometric Science of Information(GSI2017), Nov 7–9, 2017, Mines ParisTech, France.
Sei, T. (2022). Coordinate-wise transformation of probability distributions toachieve a Stein-type identity, Information Geometry, 5, 325–354.
Stein, C. (1972). A bound for the error in the normal approximation to the
distribution of a sum of dependent random variables, Proc. Sixth Berkeley
Symp. on Math. Statist. and Prob., Vol. 2, 583–602.
Villani, C. (2003). Topics in Optimal Transportation, American MathematicalSociety.
35 / 35

Mathematical Optimization and Statistical Theories using Geometric Methods

On LASSO and SLOPE estimators and theirpattern recovery
Tomasz Skalski1,21

Wrocław University of Science and Technology, Poland
2
LAREMA, University of Angers, France

Least Absolute Shrinkage and Selection Operator (LASSO) and Sorted `1 Penalized Estimator (SLOPE) are the regularization methods used for fitting
high-dimensional regression models. They allow to reduce the model dimension by nullifying some of the regression coefficients. Moreover, SLOPE allows the further reduction by equalizing some of nonzero coefficients, whichallows to identify situations where some of true regression coefficients areequal.
We shall introduce the notion of the pattern for LASSO and SLOPE and its
subdifferential-induced generalization to other convex penalized estimators,
which will be illustrated carefully in the case of the orthogonal design matrix.
This talk will present new results on the strong consistency of SLOPE estimators and on the strong consistency of pattern recovery by SLOPE when
the design matrix is orthogonal. We shall also present the relations of LASSO
and SLOPE with root system induced convex hulls.
The research was supported by a French Government Scholarship and by
Centre Henri Lebesgue, program ANR-11-LABX-0020-0.

References[1] M. Bogdan, X. Dupuis, P. Graczyk, B. Kołodziejek, T. Skalski, P. Tardivel, M. Wilczyński. Pattern Recovery by SLOPE. ArXiv 2203.12086.
[2] U. Schneider, P. Tardivel. The geometry of uniqueness, sparsity and clustering in penalized estimation. ArXiv 2004.09106.
[3] T. Skalski, P. Graczyk, B. Kołodziejek, M. Wilczyński. Pattern recovery
and signal denoising by SLOPE when the design matrix is orthogonal.
ArXiv 2202.08573.
[4] P. Tardivel, T. Skalski, P. Graczyk, U. Schneider. The Geometry ofPattern Recovery by Penalized and Structured Estimators. 2021. hal03262087.

49

50

OCAMI Reports Vol. 8 (2022)

On LASSO and SLOPE estimatorsand their pattern recovery
Tomasz Skalski
Wrocław University of Science and Technology, Poland
University of Angers, France

Osaka & on-line
2022/10/20

1 / 25

Linear regression modelLinear regression model: Y = X β + ε:

Rn : response vector
X ∈ Rn×p : design matrix
β ∈ Rp : unknown parameter vector
ε ∈ Rn : random noise term
Y ∈Noiseless case: ε = 0.
Noisy case: ε has continuous and symmetric distribution.
Goal: to estimate β.

2 / 25Mathematical Optimization and Statistical Theories using Geometric Methods

51

Ordinary Least Squares estimator

Ordinary Least Squares (Legendre, 1805, Gauss, 1809)β̂ OLS := arg minp 12 kY − Xbk22

Rb∈β̂ OLS = (X 0 X )−1 X 0 Y

Not defined when n < p.
In noisy case: with probability 1 has p pairwise different coordinates.

3 / 25Penalized estimator

Consider the following penalized estimator1β̂ := arg min kY − Xbk22 + λJ(b), where J is a norm.
b∈Rp 2β̂ is well defined when n ≥ p as well as when n < p.

The pattern of β̂ is characterized by its subdifferential ∂J .
The dual norm J ∗ is given by J ∗ (x) = sup{z 0 x : J(z) ≤ 1}.
β̂ = 0 if and only if J ∗ (X 0 Y ) ≤ 1.

4 / 25

52

OCAMI Reports Vol. 8 (2022)

Examples of penalized estimators

Ridge regression (Hoerl & Kennard, 1970)β̂ := arg minp 21 kY − Xbk22 + λ kbk2 , λ > 0Rb∈LASSO (Chen & Donoho, 1994, Tibshirani, 1996)
β̂ LASSO := arg minp 12 kY − Xbk22 + λ kbk1 , λ > 0b∈

R

SLOPE (Bogdan, van den Berg, Sabatti, Su, Candès, 2015)pP
λi |b|(i) , λ1 > 0,β̂ SLOPE := arg minp 12 kY − Xbk22 +
b∈

R

i=1

λ1 ≥ . . . , λp ≥ 0, |b|(1) ≥ . . . ≥ |b|(p)5 / 25

Least Absolute Shrinkage and Selection Operator (LASSO)

LASSO estimator (Chen & Donoho, 1994, Tibshirani, 1996) minimizes the`1 -penalized Euclidean distance between Y and Xb:
1
kY − Xbk22 +λ kbk1 ,b∈R 2β̂ LASSO := arg minp

λ > 0.

β̂ LASSO is well defined both for n ≥ p and n < p.
∂k·k1 (b) = sign(b).

6 / 25Mathematical Optimization and Statistical Theories using Geometric Methods

53

LASSO dual ball = hypercube

J ∗ (b) = kbk∞

B ∗ = B∞ (0, λ) = [−λ, λ]7 / 25Sorted `1 Penalized Estimator (SLOPE)SLOPE (Bogdan, van den Berg, Sabatti, Su, Candès, 2015) minimizes thesorted `1 penalized Euclidean distance between Y and Xb:
1β̂ SLOPE := arg minp kY − Xbk22 + JΛ (b).
b∈R 2
Sorted `1 norm: JΛ (b) :=pPi=1

λi |b|(i) , where

λ1 > 0, λ1 ≥ . . . , λp ≥ 0 and |b|(1) ≥ . . . ≥ |b|(p) .

β̂ SLOPE is well defined both for n ≥ p and for n < p.

SLOPE generalizes the previous approaches:
λ1 = . . . = λp = 0 ⇒ β̂ SLOPE = β̂ OLS ,λ1 = . . . = λp > 0 ⇒ β̂ SLOPE = β̂ LASSO .
8 / 25

54

OCAMI Reports Vol. 8 (2022)

SLOPE dual ball = Signed permutahedron P ± (Λ)
The dual of sorted `1 norm is:


|b|
|b|
+
|b|
|b|
+
·
·
·
+
|b|
(1)
(1)
(2)
(1)
(p)
JΛ∗ (b) = max
,
,...,
.
λ1
λ1 + λ2
λ1 + · · · + λp
The unit ball of JΛ∗ is the signed permutahedron P ± (Λ) :P ± (Λ) = Conv {(±λπ(1) , . . . , ±λπ(p) ) : π ∈ Sp }.

P ± (Λ) in

R2

P ± (Λ) inR3

9 / 25

Root systems and statistics

LASSO: Ap1

SLOPE: Bp

10 / 25

Mathematical Optimization and Statistical Theories using Geometric Methods

55

SLOPE pattern
Definition

Rp → Zp defined by

The SLOPE pattern is a function patt :patt(b)i = sign(bi )rank(|bi |),

i = 1, . . . , p,

where rank(|bi |) ∈ {1, 2, . . . , k} is the rank of |bi | in a set of nonzero
distinct values of {|b1 |, . . . , |bp |} (and sign(0) = 0).
Properties of patt(x):sign(patt(x)) = sign(x) (sign preservation),|xi | = |xj | =⇒ |patt(x)i | = |patt(x)j | (clusters preservation),|xi | > |xj | =⇒ |patt(x)i | > |patt(x)j | (hierarchy preservation).

Examplex = (1.2, 1.2, 5, −5, 0, 3)=⇒

patt(x) = (1, 1, 3, −3, 0, 2).
11 / 25

SLOPE vs. OLS

Theorem (Schneider & Tardivel, 2021)
For n ≥ p and ker(X ) = {0} we have:β̂ OLS − β̂ SLOPE = Proj(β̂ OLS ) on (X 0 X )−1 P ± (Λ).
For p > n:
Y − X (β̂ OLS − β̂ SLOPE ) = Proj(β̂ OLS ) on (X 0 X )−1 row(X ) ∩ P ± (Λ).

Theorem (Orthogonal design, n ≥ p)

The orthogonal projection of β̂ OLS on P ± (Λ) is equal to β̂ OLS − β̂ SLOPE .
For LASSO: proven by Ewald and Schneider (2018).

12 / 25

56

OCAMI Reports Vol. 8 (2022)

SLOPE vs. OLS
β2(-1,1)(-1,2)

(0,1)

(1,2)

(1,1)

OLS
(-2,1)

SLOPE

(2,1)

SLOPE
(-1,0)

(0,0)

OLS

(1,0)

(-2,-1)β1(2,-1)
SLOPE(-1,-1)(-1,-2)(0,-1)(1,-2)(1,-1)

OLS

Figure: β̂ SLOPE and β̂ OLS in orthogonal design: X 0 X = Ip for Λ = (2, 1)0 .

Simpler expression for SLOPE in orthogonal design: Tardivel, Servien andConcordet (2020).
13 / 25

Y (n) = X (n) β + ε(n)Consider the sequence of regression models: Y (n) = X (n) β + ε(n)with ε(n) ∼ N (0, σ 2 In ).
No assumptions on relations between ε(n) and ε(m) for n 6= m.

TheoremAssume that

lim n−1 (X (n) )0 X (n) = C > 0.
n

Let βbnSLOPE , n ≥ 1, be the SLOPE estimator corresponding to the tuning(n) (n)(n)vector Λ(n) = (λ1 , λ2 , . . . , λp )0 .
(n)

λ
a.s.
If lim 1 = 0, then βbnSLOPE −→ β.
n→∞ n
(n)
If λ0 kβk∞ > β 0 C β/2 and λ1 /n → 0,then β̂ SLOPE does not converge to β. Hence, β̂ SLOPE is not strongly
consistent for β.
14 / 25

Mathematical Optimization and Statistical Theories using Geometric Methods

57

Y (n) = X (n) β + ε(n)TheoremAssume that(n)λlim 1 = 0n→∞ n
and that there exists δ > 0 such that(n)lim inf √n→∞Then we have

λi

(n)

− λi+1n (log(n))1/2+δ

= m > 0 for i = 1, . . . , p − 1.
a.s.

patt(β̂nSLOPE ) → patt(β).

15 / 25

Application of SLOPE: signal denoising

16 / 25

58

OCAMI Reports Vol. 8 (2022)

Application of SLOPE: signal denoising

17 / 25

Application of SLOPE: signal denoising

18 / 25

Mathematical Optimization and Statistical Theories using Geometric Methods

59

Application of SLOPE: signal denoising

19 / 25

Application of SLOPE: signal denoising

20 / 25

60

OCAMI Reports Vol. 8 (2022)

Application of SLOPE: signal denoising

MSE (, ·)OLS
613.6797

LASSO-CV
426.3705

LASSO-LS
171.7957

SLOPE-LS
20.74967

Comparison of MSE between different regression methods

21 / 25

Application of SLOPE: pattern recovery

(a)

(b)

22 / 25

Mathematical Optimization and Statistical Theories using Geometric Methods

61

References
M. Bogdan, E. van den Berg, C. Sabatti, W. Su, E. J. Candès, SLOPE – Adaptive Variable
Selection Via Convex Optimization, Annals of Applied Statistics, vol.9, pp. 1103-1140,
2015.
M. Bogdan, X. Dupuis, P. Graczyk, B. Kołodziejek, T. Skalski, P. Tardivel, M. Wilczyński.
Pattern Recovery by SLOPE. ArXiv 2203.12086.
S. Chen, D. Donoho. Basis pursuit. Proceedings of 1994 28th Asilomar Conference on
Signals, Systems and Computers, 1994, pp. 41-44 vol.1.
K. Ewald, U. Schneider. Uniformly Valid Confidence Sets Based on the Lasso. EJS, vol. 12,
pp. 1358-1387, 2018.
C. F. Gauss. Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem
Ambientium, 1809.
A.-M. Legendre. Nouvelles méthodes pour la détermination des orbites des comètes, Paris,
Firmin Didot, 1805.
U. Schneider, P. Tardivel. The geometry of uniqueness, sparsity and clustering in penalized
estimation. ArXiv 2004.09106.

23 / 25

References
T. Skalski, P. Graczyk, B. Kołodziejek, M. Wilczyński. Pattern recovery and signal
denoising by SLOPE when the design matrix is orthogonal. ArXiv 2202.08573.
P. Tardivel, R. Servien, D. Concordet. Simple expression of the LASSO and SLOPE
estimators in low-dimension, Statistics 54, 340-352, 2020.
P. Tardivel, T. Skalski, P. Graczyk, U. Schneider. The Geometry of Pattern Recovery byPenalized and Structured Estimators. 2021. hal-03262087.
R. Tibshirani. Regression Shrinkage and Selection via the Lasso. Journal of the Royal
Statistical Society. Series B (Methodological), vol. 58, no. 1, [Royal Statistical Society,
Wiley], 1996, pp. 267–88.
R. J. Tibshirani, J. Taylor. The solution path of the generalized lasso. Annals of Statistics,
vol. 39, no. 3, pp. 1335-1371, 2011.
M. J. Wainwright. Sharp thresholds for high-dimensional and noisy sparsity recovery using`1 -constrained quadratic programming (Lasso). IEEE Trans. Inform. Theory, vol. 55, no. 5,
pp. 2183-2202, 2009.
P. Zhao, B. Yu. On model selection consistency of Lasso. Journal of Machine Learning
Research, vol. 7, pp. 2541-2563, 2006.
H. Zou. The Adaptive Lasso and Its Oracle Properties, Journal ASA 2006.
24 / 25

62

OCAMI Reports Vol. 8 (2022)

Domo arigato gozaimasu!

25 / 25

Appendix: Pictures from the Title Page
Meeting point of scaled B and scaled unit ball in `2 of (Y − Xb) is equalto β̂.

sign(β̂ LASSO ) = (0, +)patt(β̂ SLOPE ) = (1, 1)

26 / 25

Mathematical Optimization and Statistical Theories using Geometric Methods

63

Appendix: Subdifferential

Definition (Subgradient)
Let f :

Rp 7→ R. Then g is a subgradient of f at b if
∀h ∈ Rp f (b + h) ≥ f (b) + g 0 h.

Definition (Subdifferential)
The subdifferential ∂f (b) of f at b is the set of all subgradients of f at b.

27 / 25

Appendix: Thresholded estimator

Definition (Thresholded penalized least squares estimator)

RR

Let pen be a penalizer, X ∈ n×p , y ∈ n and λ > 0. Givenβ̂ ∈ SX ,λpen (y ), we say that û is a thresholded estimator of β̂ if∂pen (β̂) ⊂ ∂pen (û).

Definition (Thresholded LASSO)(
β̂iLASSO ,LASSO ,τβ̂i=0,

if |β̂iLASSO | > τ,otherwise.

28 / 25

64

OCAMI Reports Vol. 8 (2022)

Appendix: LASSO and SLOPE in orthogonal design
Theorem (Tibshirani, 1996)
Exact formula for LASSO in orthogonal (X 0 X = I ) design:β̂iLASSO = sign(β̂iOLS ) max{|β̂iOLS | − λ, 0}.

Theorem (Tardivel, Servien, Concordet (2020))
Let |β̂ OLS |(1) ≥ . . . ≥ |β̂ OLS |(p) . Let Ŝk :=

k
P

i=1

(|β̂ OLS |(i) − λi ). Denote apartition (k1 , k2 , . . . , ks = p) of {1, 2, . . . , p} such thatŜ −Ŝ

kk−1
ki := max{arg maxk>ki−1 { k−k}} with k0 = Ŝ0 = 0. Theni−1β̂iols · β̂islope ≥ 0 and |β̂ slope | is given by!
Ŝk1
k1 terms, . . . , (ks − ks−1 ) terms
k1+

Ŝks − Ŝks−1
ks − ks−1

! !

.

+

29 / 25

Mathematical Optimization and Statistical Theories using Geometric Methods

Likelihood Geometry of Correlation Models
Carlos Améndola
Technical University of Berlin

We present a problem where algebra appears naturally when estimatingcorrelation matrices, that is, standardized covariance matrices. Concretely,
we study the geometry of maximum likelihood estimation for correlation
matrices, which form an affine space of symmetric matrices defined by setting
the diagonal entries to one.
We study the likelihood geometry for this model and linear submodels
that encode additional symmetries. We also consider the problem of minimizing two closely related functions of the covariance matrix: the Stein’s
loss and the symmetrized Stein’s loss. Unlike the Gaussian log-likelihood,these two functions are convex and hence admit a unique positive definite
optimum.
Studying the critical points in all three settings leads to systems of nonlinear equations, and we compute some of the algebraic degree invariants
that measure the algebraic complexity of each optimization problem.
This is joint work with Piotr Zwiernik (University of Toronto, Canada).

65

66

OCAMI Reports Vol. 8 (2022)

Likelihood Geometry of Correlation Models
Carlos Enrique Améndola Cerón(Technical University of Berlin)
OCAMI: Mathematical optimization and statistical theories using geometric methods

October 20, 2022

Setup / Introduction
Sn+ real symmetric positive definite n × n matrices
Model: M ⊆ Sn+ , and Data: S ∈ Sn+

What is the ‘best’ point Σ∗ ∈ M that explains S?
Gaussian ML estimation:Σ̂ = arg max log det(Σ−1 ) − tr(Σ−1 S)Σ∈MCan be seen as minimizing the divergence I(S∣∣Σ), where−1
I(Σ1 , Σ2 ) = tr(Σ1 Σ−12 ) − log det(Σ1 Σ2 ) − n# complex critical points for generic S: ML degree
In this talk: M consists of correlation matrices, i.e. Σii = 1 ∀i
Carlos Améndola

Likelihood Geometry of Correlation Models

Mathematical Optimization and Statistical Theories using Geometric Methods

Motivating Example: Bivariate Correlations
Let M ⊂ S2+ consist of 2 × 2 correlation matrices:
1 ρ
Σ=(
)
ρ 1

K = Σ−1 =

1
1 −ρ
(
)
2
1 − ρ −ρ 1ss
S = ( 11 12 )s12 s22

where −1 < ρ < 1.
Finding the MLE corresponding to ρ̂ reduces to solving a cubic
equation [Kendall, Stuart, 1961 “Advanced Theory of Statistics”]:
ρ3 − s12 ρ2 + (s11 + s22 − 1)ρ − s12 = 0
ML degree is 3. There could potentially be three positive definite
solutions with a multimodal likelihood function `(Σ).
How often does this happen? How bad can it be?

Carlos Améndola

Likelihood Geometry of Correlation Models

A statistical perspective

4
3

2
10
50

1

2

Sample size

0

Density5

The density of the distance from the truth

0.0
Carlos Améndola

0.5

1.0

1.5

2.0

Likelihood Geometry of Correlation Models

67

68

OCAMI Reports Vol. 8 (2022)

Probability of one real critical point2

n=20.907 0.883 0.881 0.847 0.815+1.0006

0.982 0.984 0.982 0.974 0.983+0.99510+0.9900.994 0.998 0.992 0.995 0.999

14sample size

+0.9990.997 0.999 0.995 0.999 1.000+0.980+0.960+0.940+0.90018+0.8500.999 1.000 0.999 1.000 1.000
0

0.2

0.5

0.8

+0.7500.99

rho
Carlos Améndola

Likelihood Geometry of Correlation Models

Case Study: Bivariate Correlations
22
Let a = s11 +s
and b = s12 .
2Note that if S ∈ S2+ then a > 0 and ∣a∣ > ∣b∣.
d
It holds that dρ
I(S∣∣Σ) = (1−ρ2 2 )2 f (ρ), wheref (ρ) = ρ3 − bρ2 − (1 − 2a)ρ − b.
f (−1) = −2(a + b) < 0 and f (1) = 2(a − b) > 0 Ô⇒ at least one realroot in (−1, 1).
The discriminant of f is
∆f (a, b) = −4[b 4 − (a2 + 8a − 11)b 2 + (2a − 1)3 ].
f has a single real zero ⇐⇒ ∆f (a, b) < 0.
However, we are more interested in:when does f have a single critical point in (−1, 1)?
Carlos Améndola

Likelihood Geometry of Correlation Models

Mathematical Optimization and Statistical Theories using Geometric Methods

69

Likelihood Geometry for Bivariate Correlations

Carlos Améndola

Likelihood Geometry of Correlation Models

Case Study: Bivariate Correlations
Data matrix
S =(

a b
)
b a

where a > ∣b∣ > 0.
d2
2
It holds that dρ
2 I(S∣∣Σ) = (1−ρ2 )3 g (ρ), whereg (ρ) = ρ4 − 2bρ3 + 6aρ2 − 6bρ + 2a − 1.
g (−1) = 8(a + b) > 0 and g (1) = 8(a − b) > 0.
The discriminant of g is
∆g = −256(27b 6 − 27(2a2 + 6a − 5)b 4 + 9(3a4 + 36a3 − 32a2 + 8a + 1)b 2−(2a − 1)(9a2 − 2a + 1)2 ) (1)If g (ρ) > 0 for all ρ ∈ R (globally convex) Ô⇒ ∆g (a, b) ≥ 0.
However, we are more interested in:when is g nonnegative in (−1, 1)?
Carlos Améndola

Likelihood Geometry of Correlation Models

70

OCAMI Reports Vol. 8 (2022)

Convexity Analysis

Carlos Améndola

Likelihood Geometry of Correlation Models

Alternative Loss Functions
From the divergence−1
I(Σ1 , Σ2 ) = tr(Σ1 Σ−12 ) − log det(Σ1 Σ2 ) − nI(Σ1 , Σ2 ) ≥ 0 and is zero if and only if Σ1 = Σ2 .
strictly convex in Σ1 and in Σ−12
Fix S ∈ Sn+ :1

entropy loss: I(S∣∣Σ) (minimizer Σ̂ is MLE)

2

Stein’s loss: I(Σ∣∣S) (minimizer Σ̌ is dual MLE)

3symmetrized Stein’s loss:L(Σ, S) =1(I(S∣∣Σ) + I(Σ∣∣S))2(2) and (3) are strictly convex in Σ and optimizers are uniquely defined
Carlos Améndola

Likelihood Geometry of Correlation Models

Mathematical Optimization and Statistical Theories using Geometric Methods

71

Full correlation model
Let M ⊂ Sn+ be the space of all correlation matrices: Σii = 1 for all 1 ≤ i ≤ n.
First order optimality conditions give that the optimum is a correlation
matrix Σ = K −1 satisfying for each i ≠ j:1

2

entropy loss (MLE):

Kij = (KSK )ij

Stein’s loss (dual MLE):
Kij = (S −1 )ij

3symmetrized Stein’s loss:(KSK )ij = (S −1 )ij

Carlos Améndola

Likelihood Geometry of Correlation Models

Algebraic Degrees
For the bivariate correlation model n = 2, [Brownlees, Llorens-Terrazas(2020)] observed that the dual MLE can be given in closed form (solving aquadratic equation!).
From our computations, for n > 1 one has
dMLdeg(n) < MLdeg(n) < SSLdeg(n)

n

1

2

3

4

5

6

7

8

9

SSL degree

1

4

28

292

?

?

?

?

?

ML degree

1

3

15

109

1077

13695

?

?

?

dual ML degree

1

2

5

14

43

144

522

2028

8357

For n > 4, computed with the package LinearCovarianceModels.jl
how are these numbers growing?
Carlos Améndola

Likelihood Geometry of Correlation Models

72

OCAMI Reports Vol. 8 (2022)

Equicorrelation Model
The model M now consists of all Σ ∈ Sn+ such thatΣii = 1

Σij = ρ for i ≠ j.

−1
This means that ρ is restricted to n−1< ρ < 1.
−1
Let W = S . We can exploit the symmetry and set:
1
T
W =
∑ PWP
n! P∈Sn

Theorem (Am., Zwiernik (2021))
For the equicorrelation model, the dual ML degree is always 2 for everyn > 1. The dual MLE Σ̌ admits the explicit form√
1 + (n − 2)w̄ ± (nw̄ + 1)2 − 4w̄ρ̌ =.
2(n − 1)w̄where w̄ is the off-diagonal entry of W .
Carlos Améndola

Likelihood Geometry of Correlation Models

Equicorrelation Model
The model M now consists of all Σ ∈ Sn+ such thatΣii = 1

Σij = ρ for i ≠ j.

−1
This means that ρ is restricted to n−1< ρ < 1.
We can exploit the symmetry and set:
1
T
S=
∑ PSP
n! P∈Sn

Theorem (Am., Zwiernik (2021))
For the equicorrelation model, the ML degree is always 3 for every n > 1.
The MLE Σ̂ satisfies(n − 1)ρ3 + ((n − 2)(a − 1) − (n − 1)b)ρ2 + (2a − 1)ρ − b = 0.
where a, b are the diagonal and off-diagonal entries of S, respectively.
The SSL degree is always 4 for every n > 1.
Carlos Améndola

Likelihood Geometry of Correlation Models

Mathematical Optimization and Statistical Theories using Geometric Methods

Equicorrelation for n > 2

Carlos Améndola

Likelihood Geometry of Correlation Models

A statistical perspective

4
3

2
10
50

1

2

Sample size

0

Density5

The density of the distance from the truth

0.0
Carlos Améndola

0.5

1.0

1.5

2.0

Likelihood Geometry of Correlation Models

73

74

OCAMI Reports Vol. 8 (2022)

A statistical perspective

n=10

2

0.907 0.883 0.881 0.847 0.815

6

0.982 0.984 0.982 0.974 0.983

sample size

100.994 0.998 0.992 0.995 0.999

0.985 +1.000
0.981 0.950 0.871 0.854

0.997 0.999 0.995 0.999 1.000+0.999+0.995+0.995+0.990+0.990+0.980+0.9801.000 1.000 0.998 0.996 0.982
1.000 +0.960
1.000 1.000 1.000 0.999

18

18

0

0.2

0.5

0.8

0.99

0.2

Carlos Améndola

0.5

0.8

0.99

rho

rho

Likelihood Geometry of Correlation Models

SummaryRich likelihood geometry behind correlation models.
High ML degree may hint to problematic optimization, but carefulanalysis shows likelihood function is well-behaved over large regions.
Introduction of another algebraic complexity measure: SSL degree.
Dual MLE appears to behave best algebraically, how do degrees grow?
Plenty of relevant submodels (e.g. symmetries) still to be explored.
Main Reference:Améndola, C., & Zwiernik, P., Likelihood Geometry of Correlation
Models. (2021) Le Matematiche, 76(2), pp. 559 - 583.

ありがとうございました!

Carlos Améndola+0.900+0.8501.000 +0.750
1.000 1.000 1.000 1.000
0+0.960+0.9401.000 +0.900
1.000 1.000 1.000 1.000+0.8500.999 1.000 0.999 1.000 1.000+1.000+0.999+0.9401461014

sample size2

n=2

Likelihood Geometry of Correlation Models+0.750Mathematical Optimization and Statistical Theories using Geometric Methods

Mixed convex exponential families and locallyassociated graphical models
Piotr Zwiernik (University of Toronto)

Abstract
In statistical exponential families the log-likelihood forms a concave function in the canonical parameters. Therefore, any model given by convex
constraints in these canonical parameters admits a unique maximum likelihood estimator (MLE). Such models are called convex exponential families.
For models that are convex in the mean parameters (e.g. Gaussian covariancegraph models) the maximum likelihood estimation is much more complicatedand the likelihood function typically has many local optima. One solution isto replace the MLE with so called dual likelihood estimator, which is uniquelydeﬁned and asymptotically has the same distribution as the MLE. In this talkI will consider a much more general setting, where the model is given by convex constraints on some canonical parameters and convex constraints on theremaining mean parameters. We call such models mixed convex exponential
families. We propose for these models a 2-step optimization procedure which
relies on solving two convex problems. We show that the resulting estimator
has asymptotically the same distribution as the MLE. Our work was motivated by locally associated Gaussian graphical models that form a suitablerelaxation of Gaussian totally positive distributions.
(Joint work with Steﬀen Lauritzen, University of Copenhagen)

75

76

OCAMI Reports Vol. 8 (2022)

Mixed convex exponential families andlocally associated graphical models
Piotr Zwiernik
University of Toronto

This story is part of the following paper:

Lauritzen S., & Zwiernik, P., Locally associated graphical models and mixed
convex exponential families. arXiv:2008.04688.

OCAMI Meeting
21(20) October 2022

Modelling with positive
dependenceMathematical Optimization and Statistical Theories using Geometric Methods

Example: S&P 500
graphical lasso estimate of the graph:
240
349
214

121

438
356

365

20
392

26

155

157

273
281
6

426

127

217

268
173

307

178

269

36
145

51

225

67

305

144

207

287

444

52
235
338

416

37

19

139

113

447

31

449

353

230

401

441

264

119

14

254

272

111

288

250

152

8

246

433

412

35
399

245

21
397

222

7

216

62
201

34

109

64
92
105

44

196

450

90

348

345

234

87

211

96

446

135
374

60
61
123

428

425

396

68

373

327

445

312

151

343

191
10

5

104

167

23

318

197

153

377

116

147

443

162

165
141

22

125

205

375

108

40
431

328

322

212

29

166

335

341

265

279

278

82

210

193

325

133

294
110

435

236

192

319

100

137

206

274

275
136

79

391

452

38

295
163

405

129

390

160

158

83

112
395

177

297
354
232

385
350

304

103
150

347

58

370

208

384
331

321

398

323

367

164
41

430

332

258

179

255
89

283

204

262

330
260

256

180

451

360

55

437
186

300

414

386
357

248

76

174

27
380

49

394

195

420

91

130

224

361

266

342

229
302

421
159

406

46

80
344

175

18
99

39
242
43

74

393

280

241

202

253

57

388

53

50

383

409

4

311

154

97

88

24

417
314

131

427

249

86

410

142
263

271

184

376
261

182

118

161

194

339

2
231

404

146

372

72

198

424

276

310
382

25

98

183

358

346

42

65

285

73
440
422
308
333

448

363
9

189
369

117

381

324
403

1

252

413

243

340
387

371
209

317

432
259
115

56

168

101

233
32

140

439

429

316

267

289

219

315

215

226

66
45

336
419
77

218
120

134

326
213

11

221

156

293

114
290

95

334

203

402
415

301

366

149

407
378

220

94

296
75

227

143

85

3

102

257

337

247
132

223

436
106

434

389
400

138

122
188

187

200

170

309

270

48

423

282

15

169

411
359

442

244

63

30

351

12

33

379

54

126
355

16

228

306

181
277

70
368

418

299

148

190

107

329

320

185
176

352

124

47

239

291

84

286
93

303

362
13

172

69

71

128

292

251

199

59

171

298
284

17
78
408

313

28
237
238
364

81

Note: All edges green (positive partial correlations).
1

Gaussian totally positive distributions
The zero-mean Gaussian distribution
1 ?
T
fpxq “detK
expp´xK x{2qp2πqd{2Totally positive: K “ Σ´1 satisfies Kij ď 0 for all i ‰ j.
(K is an M-matrix)‚ Kij ď 0 if and only if corrpXi , Xj |XVzti,ju q ě 0.
2

77

78

OCAMI Reports Vol. 8 (2022)

A success story
In some applications it works incredibly well.
Rossell&Zwiernik describe a S&P500 dataset:‚ Our MLE gives a sparser graph and higher likelihood
than the best GLASSO estimate!
see also: Agrawal, Roy, Uhler. Covariance Matrix Estimation under
Total Positivity for Portfolio Selection, 2019.

3

However: Gene expression data
Partial correlations with negative signs additionally
penalized.
GBP

NFKBIA

MAP3K5

CEBPB

TIF
NFI

IRAK3

HIF

SLC

CXCL

MAP3K8

CEBPD

VEG

RCA
CCL2

SOCS1

NR3

TRIM8

SER
IFNGR2

ZC3

TRA
IL6

ABC

BCL3

PTGS

PTX

IL18

TNFA
IL1A

IL4

TCF

NFKBIE
PTPN2

NAM

SEL
IRF
BIR

ICA

TNI

CFL

IFIT

JAK

NFKBIZ

SOCS3

NFKB1

BDK

IFI1
IRAK1

BID

PLA2

IL6S

SPP

IL1B

MYD

CCL5

CFB

PLS

IFNGR1

C3

SIR

NFKB2

IL1RL
C1R

CAM

PLAU
LGA

STA
ISG

IL1R1

IL3

TRIM2

C1S

PTGE

REL

THB

CSF

SAM
SLA
S100A9S100A8

VDR

LCN

TRE

LY9
OLF

CXCR
IL1R2

PRK
CD1D

CTSH

BCL2

C5A

CD14

HP
CD36

TNFR

NOD

ITGAX

CEA

UCP

S100A1

SIG

CLEC7

FCG

CYT

MMP

ITGAM

PYC
TLR2

CTSS

PST
LIL

PIK
CD300L

CYBB

LYN
FPR

CD5

PTPR
NLR

CD37

SYK

CD300A

SPI
CYBA

PTPN6

TLR8

FCN
ITGB

CFP
CLEC4

LST

RAC
ARR

4

DOK

Mathematical Optimization and Statistical Theories using Geometric Methods

Motivation: Locally associated GGMs
X is associated if covpfpX q, gpX qq ě 0 for any
f, g : Rd Ñ R nondecreasing.
Pitt: A Gaussian X is associated if and only if Σ ě 0.
Gaussian graphical model: X „ Nd p0, Σq:
MpGq “ tΣ P PDd : pΣ´1 qij “ 0 for ij R Gu.
With additional positivity:
PpGq “ tΣ P PDd : Σij ě 0 for ij P Gu.
5

Estimation in laGGMs
The log-likelihood (S sample covariance matrix)
log detpΣ´1 q ´ trpSΣ´1 q
is concave in K “ Σ´1 but not in Σ.
Alternative: mixed dual estimate (MDE).
‚ MDE for mixed convex exponential families is easier toobtain and has the same asymptotics as the MLE.
6

79

80

OCAMI Reports Vol. 8 (2022)

2-stage estimation procedure
Information divergence (convex in Σ1 and in K2 ):1
1
IpΣ1 ||K2 q “ trpΣ1 K2 ´ Iq ´ log detpΣ1 K2 q.
2
2
1
2
S sample covariance, S Ñ Kp Ñ Σ̌

1. K̂ minimizer of IpS||Kq subject to K P MpGq.
2. Σ̌ minimizer of IpΣ||K̂q subject to Σ P PpGq.

Note: Σ̌ P MpGq and it is a reasonable estimator.
7

Mixed convex exponential
families

Mathematical Optimization and Statistical Theories using Geometric Methods

Regular exponential families
Exponential family E over X wrt measure νppx; θq “ exptxθ, tpxqy ´ Apθqufor θ P Θ Ď Rk ,The set of canonical parameters
"
*
ż
"
(
Θ :“ int θ P Rk :exp xθ, tpxqy νpdxq ă 8 .
X

In steep exponential families :‚ Θ convex subset of Rk ,‚ Apθq strictly convex, smooth over Θ,‚ }∇Apθq} Ñ 8 at the boundary.
8

Mixed parametrizations
The split tpxq “ pupxq, v pxqq P Rk induces splits
θ “ pθ u , θ v q P Θ,
µ “ pµu , µv q P M.
(Θ canonical parameters, M mean parameters)
Mu “ projection of M on µu
Θv “ projection of Θ on θ v
Theorem (Barndorﬀ-Nielsen, Mixed Parametrization):‚ pµu , θ v q forms an alternative parametrization.
‚ pµu , θ v q P Mu ˆ Θv (variational independence)
9

81

82

OCAMI Reports Vol. 8 (2022)

Mixed convex exponential family
Fix mixed parametrization pµu , θ v q P Mu ˆ Θv of E.
Mixed convex exponential family:‚ M1u ˆ Θ1v Ď Mu ˆ Θv
‚ M1u Ď Mu , Θ1v Ď Θv rel. closed convex subsets.
Example: Locally associated Gaussian distributions forma mixed convex exponential family.

10

Example: The Gaussian caseSuﬃcient statistics: tpxq “ ´ 12 xx T ,Canonical/mean parameters: θ “ K , µ “ ´ 12 Σ

Gradient map: ApKq “ ´ 12 log det K , ∇ApKq “ ´ 12 K ´1 .
e.g. in locally associated Gaussian graphical models:‚ Kij “ θij “ pΣ´1 qij “ 0 for ij R G , and‚ Σij “ ´2µij ě 0 for ij P G .
So this is a mixed convex exponential family.
see also Gaussian Double Markovian Distributions by Boege,
Kahle, Kretschmer, Rötger (arXiv:2107.00134)
11

Mathematical Optimization and Statistical Theories using Geometric Methods

This leads to an interesting observation:
Fix positive definite d ˆ d matrices A, B.
For any set I of indices there exists a unique positive
definite matrix Σ such that:‚ Σij “ Aij for pi, jq P I;
‚ pΣ´1 qij “ Bij for pi, jq R I.

12

Kullback-Leibler divergence
Fenchel conjugate: A˚ pµq “ suptℓpθ; µq : θ P Rk u.
Two distributions in E: one with mean parameterµp1q P M, the other with canonical parameter θ p2q P Θ.
K pµp1q , θ p2q q “ ´xµp1q , θ p2q y ` A˚ pµp1q q ` Apθ p2q q
Note: K is strictly convex both in µp1q and in θ p2q .

13

83

84

OCAMI Reports Vol. 8 (2022)

Mixed dual estimator
Mixed exponential family: pµu , θ v q P M1u ˆ Θ1v .
řn
1Suﬃcient statistics t “ n i“1 tpXpiq q “ pu, v q.

Two-step procedure:p :“ arg min K pt, θq over θ s.t. θ v P Θ1 .
(S1) θ
v
p over µ s.t. µu P M1 .
q :“ arg min K pµ, θq(S2) µu

Some properties:
q lies in the mixed convex family.
‚ Theorem: µp exists,
q exists if and only if θ‚ µ‚ if exists, it is unique (convexity),
14

Summary + bibliography + thank you!
We study submodels of exponential families where themodel constraints are convex in the mixed parameters.
Our main motivation is in local association.
The likelihood function is not concave so the MLE may
be complicated to compute.
We propose a simple and sensible alternative.
This story is part of the following paper:
Lauritzen S., & Zwiernik, P., Locally associated graphical models and mixed
convex exponential families. To appear in Annals of Statistics.
15

Mathematical Optimization and Statistical Theories using Geometric Methods

Classiﬁcation problem of invariant
q-exponential families on homogeneous spaces
Koichi Tojo
RIKEN Center for Advanced Intelligence Project

Abstract
Q-exponential family is a natural generalization of exponential family andis an important subject in the ﬁelds of information geometry and statistics.
Widely used q-exponential families such as normal distributions and Cauchy
distributions have a symmetry. More precisely, the sample space can be
regarded as a homogeneous space G/H and the family of distributions on itis G-invariant with respect to the induced G-action by pushforward. Thenthe following problem naturally arises:
Classify G-invariant q-exponential families on G/H.
I would like to talk about a strategy to solve this problem using “qdeformation” of an exponential family. Moreover, we give a new SL(2, R)invariant q-exponential family on the upper half plane.
This is a joint work with Taro Yoshino.

85

86

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Classification problem of invariant q-exponential
families on homogeneous spaces
Koichi Tojo1 , joint work with Taro Yoshino2
1 RIKEN Center for Advanced Intelligence Project, Tokyo, Japan,
2 Graduate School of Mathematical Science, The University of Tokyo

Octorber 21, 2022

1 / 36Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Contents
1

2

3

4

Introduction
ProblemMotivation
Exponential family and q-exponential family
Background
Step 1: G /H-method
Method to construct families
G -invariance of our family
Classification of G -invariant families
Step 2: q-deformation
Definition
Property
Another topic: natural projection
Questions
Example2 / 36

Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

87

ProblemExponential family and q-exponential family
Background

ProblemAim(rough)We want to know all the “good” families of distributions onimportant spaces.
Mathematically, let G be a Lie group, H a closed subgroup of G
and G /H the homogeneous space of G . Take q ∈ R.
Problem 1.1.
Classify G -invariant q-exponential families on G /H.

3 / 36Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
BackgroundA family of probability measures and machine learning
Learning by using a family of probability measures is one ofimportant methods in the field of machine learning.
Learning=to optimize the parameters in the family of
probability measures
Families of probability measures

Exponential families

4 / 36

88

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
Background

Exponential family
Exponential family
Exponential families are important subject in the field of
information geometry.
Exponential families are useful for Bayesian inference.
Exponential families include many widely used families.
Families of probability measures

Exponential families

5 / 36Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
Background

Examples (exponential families)

Table: Examples of exponential families

distributions
NormalMultivariate normal
Bernoulli
Categorical
Gamma
Inverse gamma
WishartVon Mises
Poincaré

sample sp. X

R
Rn
{±1}
{1, · · · , n}
R>0
R>0
Sym+ (n, R)
S1
H
6 / 36

Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

89

ProblemExponential family and q-exponential family
Backgroundq-exponential family
q-exponential family (q ∈ R)
q-exponential family

is a generalization of exponential family (q = 1).
is also important subject in the field of information geometry.
is useful for Tsallis statistics.
Families of probability measures
q-exponential families
Exponential families
Cauchy dists

q-Gaussian dists

7 / 36Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
Background

Relation

exponential family

q-exponential family

Amari’s
α-familyα=1

α = 2q − 1

Entropy• Shannon entropy• maximization with• Tsallis entropy• extremization withexpected value constraint

q-expected value constraint

8 / 36

90

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
BackgroundDefinition of q-exponential family
X : manifold, R(X ): the set of all Radon measures on X .
Definition 1.2 (q-exponential family).
P ⊂ R(X ) is an q-exponential family on X if there exists a triple(µ, V , T ) such that1
2
3
4µ ∈ R(X ),V is a finite dimensional vector space over R,
T : X → V , x 7→ T (x) is a continuous map,

For any p ∈ P, there exists θ ∈ V ∨ such that

dp(x) = cθ−1 expq (−⟨θ, T (x)⟩)dµ(x),∫where cθ = x∈X expq (−⟨θ, T (x)⟩)dµ(x) (normalizing
constant).

We call the triple (µ, V , T ) a realization of P.
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

9 / 36

ProblemExponential family and q-exponential family
BackgroundDefinition of expq
For q ∈ R, we put Iq := {x ∈ R | (1 − q)x + 1 > 0}.
Definition 1.3.
The map expq : Iq → R>0 is defined by
expq x :={ex((1 − q)x + 1)

(q = 1),1
1−q(q ̸= 1).

Remark 1.4.
expq is defined as the inverse map of the q-logarithm function
lnq : R>0 → R
{
∫ x
ln x
(q = 1)1
lnq x :=dt=1q1−q − 1) (q ̸= 1).
1 t
1−q (x
10 / 36Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

91

ProblemExponential family and q-exponential family
Background

Graph of expq
expq x :=e x(q = 1),1((1 − q)x + 1) 1−q(q ̸= 1).

q<0

q=0

q=1

0<q<1

q>1

11 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
Background

Example: a family of normal distributions
Example 1.5.
The following family of normal distributions is an exponential
family on R (q = 1):
{
(
) }
1
(x − m)2
P := √
exp −
dx
2σ 2
2πσ 2
(σ,m)∈R>0 ×R
1µ =Lebesgue measure,2

V = R2 ,

3

( 2)
x
.
T : X = R → R2 , x 7→
x(µ, V , T ) is a realization of P.
12 / 36

92

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
Background

Example: a family of Cauchy distributions
Example 1.6.
The following family of Cauchy distributions is a 2-exponential
family on R:
{
}
1
γ
P :=
π (x − x0 )2 + γ 2 (γ,x0 )∈R>0 ×R
1µ =Lebesgue measure,2

V = R2 ,3

T : X = R → R2 , x 7→

( 2)
x
.
x

(µ, V , T ) is a realization of P.
13 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
Background

Background
Remark 1.7.
By definition, there are too many q-exponential families.
Only a small part of them are widely used.

q-Exponential families enlarge
Normal dists
Gamma dists
Normal dists
Gamma distsvon Mises dists

Bernoulli dists

q-Gaussian distsvon Mises dists
Bernoulli dists

Widely used exp. families

q-Gaussian dists

14 / 36Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

93

ProblemExponential family and q-exponential family
Background

Motivation
We can expect there exist “good” q-exponential families.
We want a framework to understand “good” q-exponential families
systematically.

q-exponential families
good↓Gauss dists
Gamma distsvon Mises dists
Bernoulli dists

Cauchy dists

Widely used q-exp. families
15 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
Background

Observation 1.8.
Useful q-exp. families have the same symmetry as the sample
spaces.
Sample space : homogeneous space G /HFamily : invariant under the induced G -action

normal dist

R× ⋉ R-action(scaling
andtranslation)

normal dist−→

G /H = (R× ⋉ R)/R× ≃ R
16 / 36

94

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

ProblemExponential family and q-exponential family
Background

Strategy

Problem 1.1 (again)
Classify G -invariant q-exponential families on G /H.
Step 1 Classify G -invariant exponential families on G /H by using
G /H-method.
Step 2 Classify G -invariant q-exponential families on G /H by
q-deformation of G -invariant exponential families on G /H.

17 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Method to construct families
G -invariance of our family
Classification of G -invariant families

G /H-method
We proposed a method to construct exponential families.
The method generate many well-known families.
Families obtained by the method can be classified.

Exponential families

Our method
Normal dists
Gamma distsvon Mises dists
Bernoulli dists

Poisson dist

Widely used exp. families
18 / 36Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

95

Method to construct families
G -invariance of our family
Classification of G -invariant families

G /H-method: overview

G /H-method = a method to construct a family of
probability measures on G /H froma finite dim. real representation ρ : G → GL(V ),a nonzero H-fixed vector v0 ∈ V .
See [TY18, TY19, TY20] for the details.

19 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Method to construct families
G -invariance of our family
Classification of G -invariant families

Examples obtained by our method
Table: Examples and inputs (G , H, V , v0 ) for them

distributions sample sp. X
G
H
V
v0
×
×
NormalR
R ⋉R
R
Sym(2, R)E22nnMulti. normalR
GL(n, R) ⋉ R GL(n, R)Sym(n + 1, R)En+1,n+1
Bernoulli
{±1}
{±1}
{1}
Rsgn
1
Categorical {1, · · · , n}
Sn
Sn−1
W
w
Gamma
R>0
R>0
{1}
R
1
Inverse gamma
R>0
R>0
{1}
R−1
1
+
WishartSym (n, R) GL(n, R) O(n) Sym(n, R)
In
1
2
Von Mises
S
SO(2)
{I2 }
R
e1
Poincaré
H
SL(2, R) SO(2)I2∑n Sym(2, R)nHere W = {(x1 , · · · , xn ) ∈ R | i=1 xi = 0},w = (−(n − 1), 1, · · · , 1) ∈ W .

20 / 36

96

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Method to construct families
G -invariance of our family
Classification of G -invariant families

Example: Poincaré dists on the upper half plane
Upper half plane H := {z = x + iy ∈ C | y > 0} admits the linear
fractional transformation of SL(2, R).
⇝G = SL(2, R), H = SO(2), X := G /H ≃ H.
geodesics

Low dimensional representation:
ρ : SL(2, R) → GL(Sym(2, R)),ρ(g )S = gS tg (S ∈ Sym(2, R)).
v0 := I2 .
(
}
{ 2D
)
De
a(x 2 + y 2 ) + 2bx + c dxdy (
)
exp −
⇝
a b
π
y
y2
∈Sym+ (2,R)
b c√Here D = ac − b 2 .
Higher dimensional cases:
We obtain new families by G /H-method.

21 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Method to construct families
G -invariance of our family
Classification of G -invariant families

P is a G -invariant exponential family
Theorem 2.1.
Any family obtained by our method is a G -invariant exponential
family on G /H.

normal dist

R× ⋉ R-action(scaling
andtranslation)

normal dist−→

We obtain a family with the symmetry of G /H !
22 / 36Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

97

Method to construct families
G -invariance of our family
Classification of G -invariant families

Question
Conversely,
Question 2.2.
Are any G -invariant exponential families on G /H obtained by our
method?
⇝ Yes, under a mild assumption.
⇝ Roughly speaking,{G -invariant exponential family on G /H}“=”{family on G /H obtained by G /H-method}

23 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Method to construct families
G -invariance of our family
Classification of G -invariant families

Answer to the question
Setting 2.3.
P := {pθ }θ∈Θ is a G -invariant exponential family on G /H. Here Θis the parameter space.
Theorem 2.4.
Assume1

G /H admits a nonzero relatively G -invariant measure,2

Θ is open.

Then, P is a subfamily of a certain family obtained by
G /H-method.
For the details, see our paper [TY20].
24 / 36

98

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Method to construct families
G -invariance of our family
Classification of G -invariant families

Classification of G -invariant exponential families

Let us consider an important homogeneous space G /H such as asphere and a hyperbolic space, more generally symmetric spaces.
Step 1
Classify G -invariant exponential families on G /H.
By Theorem 2.4, this problem above is reduced to the following:
Question 2.5.
Classify families obtained by G /H-method on G /H.

25 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Definition
Property

q-deformation of exponential family
q-deformation is a method to construct a q-exponential family
from an exponential family with its realization.
Definition 3.1.
Let P be an exponential family on X and (µ, V , T ) realization of
P. Putd q̃θ (x) := expq (−⟨θ, T (x)⟩)dµ(x)Θ := {θ ∈ V∨

qθ := cθ−1 q̃θ ,Pq := {qθ }θ∈Θ(θ ∈ V ∨ , x ∈ X )| − ⟨θ, T (x)⟩ ∈ Iq for any x ∈ X ,∫cθ :=d q̃θ (θ ∈ Θ)∫X

d q̃θ < ∞}X

Then Pq is a q-exponential family on X . We call Pq a
q-deformation of exponential family (P, (µ, V , T )).

26 / 36Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

99

Definition
Property

ExampleExample 3.2.
The family of Cauchy distributions is obtained by 2-deformation of
the family of normal distributions.
P: family of normal distributionµ: Lebesgue measure,
V = R2 ,T : X = R → R2 , x 7→

( 2)
x
.
x

⇝ P2 : the family of Cauchy distributions.
27 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Definition
Property

Property of q-deformation

Let X := G /H be a homogeneous space admitting nonzero
relatively G -invariant measure and q ∈ R.
Proposition 3.3.
Let P be a G -invariant exponential family on X . Then, there
exists a realization (µ, V , T ) of P such that µ is a relatively
G -invariant measure on X . Moreover, If q > 1 and P is full , thenthe q-deformation Pq of (P, (µ, V , T )) is G -invariant
q-exponential family on X .

28 / 36

100

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Definition
Property

Question
Conversely,
Question 3.4.
Are any G -invariant q-exponential families on G /H obtained by
q-deformation of some exponential family?
⇝ Yes if q > 1 under a mild assumption.
⇝ Roughly speaking,{G -invariant q-exponential family on G /H}“=”{q-deformation of G -invariant exponential family on G /H}

29 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Definition
Property

Answer to the question
Setting 3.5.
Pq = {pθ }θ∈Θ is a G -invariant q-exponential family on G /H
(q > 1).
Theorem 3.6.
Assume1

G /H admits a nonzero relatively G -invariant measure,2

Θ is open.

Then, Pq is a subfamily of a q-deformation of a certain G -invariant
exponential family with a relatively G -invariant base measure.

30 / 36Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

101

Definition
Property

Classification of G -invariant q-exponential families
G /H: important space
Step 2
Classify G -invariant q-exponential families on G /H by using
q-deformation.
By Theorem 3.6, this problem above is reduce to the following:
Question 3.7.
Classify families obtained by q-deformation of G -invariant
exponential families.

31 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Definition
Property

New family of distributions on the upper half plane
Theorem 3.8 (q-deformation of the family of Poincaré
distributions).
Let q ∈ [1, 2). The following family of distributions is
SL(2, R)-invariant q-exponential family on the upper half plane.
{()}
2 + y 2 ) + 2bx + c
dxdy
a(x


cθ−1 expq −yy 2 θ:=a b ∈Sym+ (2,R)
b c√π(expq (−2D))2−qcθ :=, D := ac − b 2 .
(2 − q)D

32 / 36

102

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Questions
ExampleFamily on a group and the natural projection

G : Lie group, H: closed subgroup of G ,π : G → G /H, g 7→ gH natural projectionπ∗ : P(G ) → P(G /H) pushforward
Question 4.1.
What kind of families can we obtain by the pushforward π∗ of
G -invariant exponential family on G ?
Is pushforward of exponential family also exponential family?
Is pushforward of G -invariant family also G -invariant?

33 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Questions
ExampleFamily on a group and the natural projection
G : Lie group, H: closed subgroup of G ,π : G → G /H, g 7→ gH natural projectionπ∗ : P(G ) → P(G /H) pushforward
Question 4.1.
What kind of families can we obtain by the pushforward π∗ of
G -invariant exponential family on G ?
Is pushforward of exponential family also exponential family?
⇝No, in general.
Is pushforward of G -invariant family also G -invariant?
⇝Yes.

33 / 36Mathematical Optimization and Statistical Theories using Geometric Methods
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

103

Questions
ExampleExample: exp. family on aﬃne transformation group
G := R>0 ⋉ R, α ∈ R,1ρα : G → GL(3, R), ρα (a, b) := bv0 := t(1, 0, 0).

b2
2


0 0
1 0 diag(aα , aα+1 , aα+2 ),b 1

Proposition 4.2.
If α ̸= 0, by applying G /H-method to (ρα , v0 ), we get a
R>0 ⋉ R-invariant exponential family on R>0 ⋉ R as follows:
{
}
)
√ (
|α| u det D rα
2α(r + 21 )−1√exp(−a (s + tb + ub ))adadbuπΓ(r )(r ,S)∈Θ
Here, D :=

(

s
t
2

t
2

u

)
, (a, b) ∈ G and Θ := R>0 × Sym+ (2, R).
34 / 36

Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Questions
ExamplePushforward of the obtained family
P :=

{

}
√
|α| u ( det D )rα
2α(r + 12 )−1√exp(−a
(s
+
tb
+
ub
))a
dadb
u
πΓ(r )(r ,S)∈Θ

G := R>0 ⋉ R, H := R>0 , π : G → G /H ≃ R, π∗ : P(G ) → P(R).
Proposition 4.3.
The family π∗ P on R is given as follows:
{
√
(
)}
1
2
Γ( q−1)q−11(b − m)expq −
√
3−q
2 σ
2σ 2
πΓ( 2(q−1) )(q,m,σ)∈(1,3)×R×R>0

Remark 4.4.
Each distribution is a q-Gaussian distribution.
The family does not depend on α ∈ R× .
35 / 36

104

OCAMI Reports Vol. 8 (2022)
Introduction
Step 1: G /H-method
Step 2: q-deformation
Another topic: natural projection

Questions
ExampleReferences[TY18] K. Tojo, T. Yoshino, A method to construct exponential
families by representation theory, arXiv:1811.01394v4, toappear in Info. Geom.
[TY19] K. Tojo, T. Yoshino, On a method to construct exponential
families by representation theory, GSI2019, Lecture Notes in
Computer Science, vol 11712, 147–156 (2019).
[TY20] K. Tojo, T. Yoshino, Harmonic exponential families on
homogeneous spaces, Info. Geo. (2020).
https://doi.org/10.1007/s41884-020-00033-3

36 / 36Mathematical Optimization and Statistical Theories using Geometric Methods

105

Adaptive shrinkage of singular values for a low-rank matrix mean when a covariancematrix is unknown
Yoshihiko Konno
Department of Mathematics, Osaka Metropolitan University
Assume that
" m,#n, p are positive integers" such# that" min{m,# n} ≥ p and that we observeXXΞ
a matrix
which is modeled as=+ E where Ξ is an m × p
YY
0n×p
non-random matrix(unknown and its rank may be less than min{p, m}), E is an (m +n) × p error matrix(unobservable) whose rows are identically distributed as Np (0p , Σ), a
p-variate real normal distribution with zero mean vector and covariance matrix Σ. We
assume that Σ is a p × p positive-definite and unknown matrix.
We consider the problem of estimating Ξ under a low-rank mean matrix condition, i.e.,rank Ξ = r < p;r is unknown©
ªb Ξ| Σ) = tr (Ξb − Ξ)> (Ξb − Ξ)Σ−1 , where Ξb := Ξ(X,bunder a loss function L(Ξ,Y ) is>an estimator of Ξ. Here A and trA stand for the transpose and the trace of a squareb Ξ| Σ) is given by the expected value of the loss
matrix A. The risk function of R(Ξ,function where the expectation is taken with respect to the joint distribution of (X, Y ).
We give Steins’s unbiased risk estimate for estimators of the form
b=
Ξ

µX
p
j=1

hj (`j )uj v >j¶(Y > Y )1/2 .

Here hj : [0, ∞) → [0, ∞), (j = 1, 2, . . . , p) are absolutely continuous functions and
U LV > is the singular value decomposition of X(Y > Y )−1/2 where U = (u1 , u2 , . . . , up )is an m × p matrix such that U > U = I p (the p × p identity matrix), V = (v 1 , v 2 , . . . , v p )is a p × p orthogonal matrix, and L is a p × p diagonal matrix whose j-th diagonal
element is given by `j . Note that we may assume that `1 > `2 > · · · > `p > 0 (almost
everywhere) with out loss of generality. Based on SURE formula, we propose an adaptive
soft-theshholding rule to the singular values `1 , `2 , . . . , `p . Furthermore, the results above
are extended to the complex normal distribution setup.

106

OCAMI Reports Vol. 8 (2022)

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Adaptive shrinkage of singular values of a low-rankmean matrix when a covariance matrix is unknown
Yoshihiko KONNO
Osaka Metropolitan Univeristy/JWU

Workshop
Mathematical optimazaiton and statistical theories using geometric methods

20-21 Octorber 2022

1 / 31MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

1

MANOVA model and its canonical mode2

Problem set-up

3

Mean matrix estimation when a covariance is known

4

Mean matrix estimation when a covarianc matrix is unknown

5

Concluding remarks

2 / 31

Mathematical Optimization and Statistical Theories using Geometric Methods

107

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

1

2

3
4

1

The reconstruction of a low-rank matrix from its noisy
observation is useful in many applications. This problem isreformulated into a constrained nuclear norm minimization
problem (regularized problem).
An important ingrident of this problem is how to choose aregularization parameter based on data. Usually the data is
independently and identically distributed with unknown
variance.
(1) The discrepacy principle approach, (2) Stein’s Unbiased
risk estimator(SURE) approach.
Inspired by approach(2) we consider the problem of
estimating a low-rank matrix mean in MANOVA(Mulitivariate
Analysis of Variance) setting 1 when a positive-definite
covariance matrix of error is unknown.

We have data for unknown covariance matrix. The distribution of this data ismean-zero.

3 / 31MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

MANOVA model and its canonical model
Let m, n, p ∈ N such that min(m, p ) ≥ p . Consider amultivariate regression modelW= |{z}A
Err ,B + |{z}|{z}|{z}(m+n)×p(m+n)×m m×p(m+n)×pwhere A is a known design matrix of full rank, B is an unknown
regression matrix of rank r (< min(m, p ) and r is unknown), and
Err is an unobservable error matrix. Here rows of Err are
independently and identically distributed as Np (0p , Σ) where Σ isa p × p positive-definite unknown matrix.

4 / 31

108

OCAMI Reports Vol. 8 (2022)

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Notation1 e ⊤1 e ⊤
2
Err = 
..

.

⊤
e





 : (m + n) × p ,




m +n

2





vec(Err) := 



e1
e2
..
.
e m +nwhere e j ’s are independently and identically distributed as
Np (0p , Σ) (j = 1, 2, . . . , (m + n)).
Write









[{
}{
}⊤ ]
COV(Err) = E vec(Err − E[Err] vec(Err − E[Err]= I m+n ⊗ Σ,

Err ∼ N(m+n)×p (0(m+n)×p , I m+n ⊗ Σ).

5 / 31MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

1Let

P = (A ⊤ A )−1/2 A ⊤ : m × (m + n)

and take P ⊥ : n × (m + n) s.t.

P ( P ⊥ ) ⊤ = 0 m ×n

and P ⊥ (P ⊥ )⊤ = I n ,Note that[2(

Put Ξ := A ⊤ A

[

X
Y

]

:=

[

P
P⊥

)1/2

P
P⊥

]

]

[P ⊤ , (P ⊥ )⊤ ] = I m+n .

B andW ∼ N(m+n)×p

([

Ξ
0n×p

]

)
, I m +n ⊗ Σ .

6 / 31

Mathematical Optimization and Statistical Theories using Geometric Methods

109

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Problem set-up
Assume that min{m, n} ≥ p and that

m
n

(

p
X
Y

)

=

(

)
Ξ
+ E;
0

[

p


 ←→  ←→  . E =  ..  ←→ 


←→

]
X
is observation and Ξ is an m × p non-random
whereY
matrix(unknown ) of rank r < p , E is an (m + n) × p error
matrix(unobservable) whose rows are identically distributed as
Np (0, Σ). Here Σ is a p × p positive-definite and unknown matrix.

7 / 31MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

We consider the problem of estimating Ξ under a low-rank mean
matrix condition, i.e.,rank Ξ = r < p ;r is unknown

under a loss fucntion and its riskb Ξ) = tr {(Ξb − Ξ)Σ−1 (Ξb − Ξ)⊤ } =: ∥Ξb − Ξ∥2LΣ (Ξ,F, Σandb Ξ) = E[LΣ (Ξ,b Ξ)]RΣ (Ξ,b is an estimator based on (X , S ). Here
where ΞS = Y ⊤ Y ∼ Wp (Σ, n), which is the Wishart distribution with thedegree of freedom n and the scale matrix Σ.

8 / 31

110

OCAMI Reports Vol. 8 (2022)

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Mean matrix estimation when a covariance is known
Assume that m, p are postive integers s.t. m ≥ p .
Let

 ⊤ 
 z 1  ⊤  z 
Z =  .2 
 .. 


T 
zm

be an m × p data matrix whose row vectors are independently
distributed asz i : p × 1 ∼ N(eξi , σ2 I p ),(i = 1, 2, . . . , m)⊤e := (eHere Ξ
ξ1 , . . . , eξm ) is unknown but σ > 0 are known.

9 / 31MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

We assume that low-rank mean matrix condition, i.e.,( )e = r < p;rank Ξ

r is unknown.

e under a loss fucntion
Consider the problem of estimating Ξand its risk
and

b

b

b

b

e Ξ)
e = tr {(Ξe − Ξ)(e Ξ
e − Ξ)e ⊤ } =: ∥Ξe − Ξ∥e 2
L1 (Ξ,
F

b

b

b

e Ξ)
e = E[L1 (Ξ,e Ξ)].
e
R1 (Ξ,

e is an estimator based on Z .
Here Ξ
tr A and A ⊤ stand for the traace and the transpose of amatrix A , respectively.

√ ()∥A ∥F := tr A ⊤ A , the Frobenius norm of a matrix A .

10 / 31

Mathematical Optimization and Statistical Theories using Geometric Methods

111

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Eckart-Young approximation theorem
Singular Value Decomposition: We can assume that m ≥ p
without loss of generality. Decompose Z asZ

= ULV ⊤ ;L

= diag(ℓ1 , ℓ2 , . . . , ℓp ) with ℓ1 ≥ ℓ2 ≥ · · · ≥ ℓp ≥ 0

U = (u1 , . . . , up ), V = (v 1 , . . . , v p )where ui ∈ Rm , v i ∈ Rp (i = 1, . . . , p ) s.t.
U ⊤ |{z}U = V ⊤V = Ip .
m ×p

The total least squres (TLS) pseudo estimator is given bybe TLS =
Ξr∑

ℓi ui v ⊤
.
i

⇐⇒

b
e TLS ∈
Ξe − Z ∥2 .
argmin∥Ξ( )F

e rank Ξe ≤rΞ:(i =1)
Notaton σj A > 0 (j = 1, 2, . . . , r ) are non-zero singularvalues of a matrix A with r = rank (A ).

11 / 31

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Regularization approach
We consider an estimator which minimizes the penalized leastsquares criterionMat(m, p ; R) ∋ Ξ 7→
( )

1
2

( )
∥Z − Ξ∥2F + penλ Ξ ∈ [0, ∞)where penλ · (≥ 0) is a penality function of Ξ and λ (≥ 0) is
a tuning parameter.

12 / 31

112

OCAMI Reports Vol. 8 (2022)

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Examples of penalities: For a positive λ > 0,( )( )⋆ penλ Ξ = λrank Ξ=⇒ a hard-theshholding rule, i.e.,SVHTλ (Z ) =where 1l{event} =

{

j =1

if event is true,.
otherwise

1
0

p
∑
( )
( )
σΞ
⋆ penλ Ξ = λ∥Ξ∥1 := λ
j =1

{

( )

where (σ Ξ , uj , v j )jp∑}min(m, p )j =1

jℓj 1l{ℓj ≥ λ}uj v ⊤
,
j

( ( )
)
σ Ξ : SV’s of Ξjis a system of singular values of Ξ=⇒ a soft-thresholding rule, i.e.,p∑()ℓj − λ 1l{ℓj ≥ λ}uj v ⊤SVSTλ (Z ) =.
j

13 / 31

j =1

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

A hard-shreshholding rule
Assume that σ2 is known.
Solve

]

[1

SVHTλ (Z ) = argmin ∥Ξ − Z ∥2F + λ rank (Ξ)
2Ξwhere λ > 0 is a tuning scalar parameter.
Then the solution is given by
SVHTλ (Z ) =

p
∑
j =1

ℓj 1l{ℓj ≥ λ}uj v ⊤; 1l{ℓj ≥ λ} =
j

{

1
0

4 √
The optimal shreshholding is √p σ when p = m.
3
(See Donoho and Garvish (2017, IEEE, Trans. Inform
Theory).

ℓj ≥ λotherwise

14 / 31

Mathematical Optimization and Statistical Theories using Geometric Methods

113

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Steps to obtain an adaptive thresholding esitmator
1

2

Solve regularized minimizaton problem

b
eλ ∈
Ξ

{
( )}
2
arg min∥Z − Ξ∥F + penλ Ξ . .

Ξ∈Mat(m ,p ; R)

(

b )

eλ :
Calculate SURE if possible a closed form of Ξ[(b)( b )]eλ , Ξ
e = E SURE Ξ
eλ
R1 Ξ

3

(b )

e λ is a function of λ and observable data.
Note that SURE Ξ
Solve minimization problem(b )be λ =⇒ Ξeb.
b
λ ∈ SURE Ξ
λ15 / 31

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Remarks
1

This method works for the soft-thresholding rule. See Cándes
el al. (2013).

2

SURE does not work for the hard-thresholding rule since
Stein’s identity, integration-by-parts formula with respect tomultivariate normal distribution, fails for the hard-thresholding
rule becuase of discontinuity of estimator.

16 / 31

114

OCAMI Reports Vol. 8 (2022)

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

A soft-thersholding ruleCèndes et al. define an adaptive soft-shreshholding rule
based on SURE:
p
p
∑
∑
⊤
SVSTλ (Z ) =(ℓi − λ) 1l{ℓi ≥ λ}uj v j =:(ℓj − λ)+ uj v ⊤
j
j =1

j =1

(1)

which is obtained from

{
}
p
∑
1
λj
min ∥Z − Y ∥2F + λ
Y 2
j =1

Y = SVSTλ (Z ).

The parameter λ in (1) is selected by minimizieng SURE,
Stein’s unbiased risk estimate for (1).

17 / 31
MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Gaussian integration-by-parts (=Stein’s identity) and a bit ofalgebraic calculation lead to
R1 (SVSTλ , Ξ) = E[SURE(SVSTλ )(Z )],2

SURE(SVSTλ )(Z ) = −mp σ +

p
∑
j =1

min{ℓj2 , λ2 }+2σ2 div(SVSTλ (X )),
)
p (
p
∑
∑
λ
div(SVSTλ (Z )) = (m − p )1−+1l{ℓj > λ}
ℓ
j
+
j =1
j =1
+2

p ∑
p
∑
ℓj (ℓj − λ)+j =1 k , iwhenever ℓ1 > ℓ2 > · · · > ℓp ≥ 0.

ℓ2 − ℓ2jk

18 / 31

Mathematical Optimization and Statistical Theories using Geometric Methods

115

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

An adaptive estimator is given byp∑
SVSTbλ (Z ) =(ℓi − bλ)+ uj v ⊤
,
i

(2)

j =1

b
λ ∈

 p

∑

arg min min{ℓi2 , λ2 } + 2σ2 div(SVSTλ (Z )) .
λ≥0

i =1

Numerical evaluation of the risk of (2) was carried out by
Candés et. al.

e is close to
But it is not clear if R1 (SVSTbλ (Z ), Ξ)( )bTLS (Z ), Ξ)e for ∀Ξe s.t. rank Ξe ≤ r < min(m, p ).
R1 (Ξ19 / 31
MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Mean matrxi estimation when a covarianc matrix is unknown
Assume that min{m, n} ≥ p and that

m
n

(

p
X
Y

)

=

(

)
Ξ
+ E;
0

p


 ←→  ←→  .. 
E =  .  ←→ 
←→The m × p mean matrix Ξ is of rank r < p
The error E is an (m + n) × p error matrix(unobservable)
whose rows are identically distributed as Np (0, Σ).
The covariance matrix Σ is a p × p positive-definite andunknown.

20 / 31

116

OCAMI Reports Vol. 8 (2022)

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

We consider the problem of estimating Ξ under low-rankmean matrix condition, i.e.,rank Ξ = r < min(m, p );r is unknown.

A loss fucntion and its risk are given byb Ξ) = tr {(Ξb − Ξ)Σ−1 (Ξb − Ξ)⊤ } =: ∥Ξb − Ξ∥2LΣ (Ξ,F, Σandb Ξ) = E[LΣ (Ξ,b Ξ)]RΣ (Ξ,b is an estimator based on (X , S ).
where ΞS = Y ⊤ Y ∼ Wp (Σ, n), which is the Wishart distribution withthe degree of freedom n and the scale matrix Σ.

21 / 31
MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

To derive a class of estimators, first assume that Σ is known.
Then we havee I m ⊗ I p ),X Σ−1/2 ∼ Nm×p (Ξ,e given bywhich leads to an estimator of Ξ
b
e TLS
Ξ

∈ arg min ∥X Σ−1/2 − Ξ∥F 2
rank Ξ≤re = ΞΣ−1/2
Ξ

=⇒

b
b=Ξ
e TLS Σ1/2 .
ΞHence we consider a class of estimators of the form

 p


∑
⊤
 S 1/2 ; X S −1/2 = ULV ⊤bH = 
Ξ
h
(ℓ
)
u
v
i
i
i


i 
i =1

where L = diag(ℓ1 , . . . , ℓp ), H = diag(h1 , . . . , hp ),U = (u1 , . . . , up ) and V = (v 1 , . . . , v p ) s.t.
U⊤U = V ⊤V = Ip .

22 / 31

Mathematical Optimization and Statistical Theories using Geometric Methods

117

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Regularized minimization problem
Known Σ case: For λ ≥ 0,

Mat(m, p ; R) ∋ ΞΣ−1/2
e −1/2 ∥F 2 + 2λ∥ΞΣe −1/2 ∥1
7→ ∥X Σ−1/2 − ΞΣbe of aUnknonw Σ case: For λ ≥ 0, find a minimizer Ξ
regularized minimization problem

andeMat(m, p ; R) ∋ Ξ

e F 2 + 2λ∥Ξ∥e 1
7→ ∥X S −1/2 − Ξ∥bb=Ξe S 1/2 =
Ξ

(∑
j =1

)
ℓj (ℓj − λ)+ uj v ⊤S 1/2
jwhere {(ℓj , uj , v j )} is a system of singular values of X S −1/2 .

23 / 31

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarksIfc(j = 1, 2, . . . , p );ℓj
c is a known positive constant,hj (ℓj ) = ℓj −then it results in the Efron-Morris estimator which is given by[]bH = X S −1/2 I p − c {(X S −1/2 )⊤ (X S −1/2 )}−1 S 1/2
Ξ= X − cX {X ⊤ X }−1 S .

On the other hand, Tsukuma and Kubokawa (2015)
considered estimators of the formbT = X − UT U ⊤ X
Ξwhere T = diag(t1 (ℓ2 ), . . . , tp (ℓp2 )) and X S −1/2 = ULV ⊤1with m × min(m, p ) matrix U s.t. U ⊤ U = I min(m, p ) .

24 / 31

118

OCAMI Reports Vol. 8 (2022)

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Recall thatX S −1/2 = ULV ⊤ ⇐⇒ L −1 U ⊤ X = V ⊤ S 1/2 .
From a simple calculation we get

bH = UHV ⊤ S 1/2 = UHL −1 U ⊤ X = UL −1 HX .
Ξ√
If we set I p − T = L −1 H (tj (x ) = hi ( x )), then we have

From this we can see thatbT .
bH = Ξ
ΞbH , Ξ) = LΣ (ΞbT , Ξ).
LΣ (Ξ25 / 31

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Furthermore, using the result due to Tsukuma andKubokawa (2015), we havebT , Ξ) = E[SURE(T )];
RΣ (Ξ

p [
∑
SURE(T ) =m + a ℓj2 tj2 − 2bti − 4ℓj2 tje
tj − 4ℓi2e
tj
j =1

4 2

4 2

2

p ℓ t −ℓ t
∑
j j
k k

ti
a, b

2

p ℓ tj − ℓ tk ]
∑
j
k
−2
−4
;
2
2
2
2
ℓ
−
ℓ
ℓ
−
ℓ
k ,j
k ,j
j
k
j
k


hi (ℓj )h (ℓj ) 
1 e′′ ,h (ℓj ) += 1−;
ti = −
j
2
ℓj
ℓj
2ℓ
i

:

known positive constants.

26 / 31

Mathematical Optimization and Statistical Theories using Geometric Methods

119

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Then we have an adaptive soft-thresholding rule

bb
Ξ
λ

 p

∑


−1/2
1/2
⊤
b
= : SVSTbλ (X S)S=  (ℓj − λ)+ uj v j  S 1/2
j =1

where bλ = argminλ≥0 SURE(SVSTλ )(X S −1/2 );SURE(SVSTλ )(X S

−1/2

)=

p [
∑
j =1

m + a ℓj2 tj2 − 2bti − 4ℓj2 tje
tj

4 2− 4ℓj2etj − 2

tj = 1 −4 2

p ℓ t −ℓ t
∑
j j
k k
k ,j

(ℓj − λ)+ℓ2 − ℓ2jk∑ ℓi2 ti − ℓi2 ti ]
−4
;
2
2
ℓ
−
ℓ
j ,i
i

j

(j = 1, . . . , p );ℓj(ℓj − λ)+ 

−
2
 .
e
tj = −(2ℓj ) 1l{ℓj > λ} +ℓj27 / 31

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Special case

Σ = σ2 I p where σ is postive but unknown.
Let s 2 = tr(Y ⊤ Y )/p .
Then an adaptive
soft-thresholding
rule for this case is given
)
∑p (
⊤
2
b
b
by Ξbλ =ℓi − λs + uj v ; , X = ULV ⊤ , withj =1
jb
λ = argminSURE(SVSTλ )(X ) andλ≥0

SURE(SVSTλ )(X ) =

p [
∑
j =1

4 2

ms 2 + a ℓj2 tj2 − 4ℓje
tj − 2(+s a ℓj2 tj2 − 4ℓj2 tje
tj − 42

4 2

p ℓ t −ℓ t
∑
j j
k k
k ,j

ℓ2 − ℓ2jk2
p ℓ tj − ℓk tk )]
∑
j
k ,j

ℓ2 − ℓ2jk

28 / 31

120

OCAMI Reports Vol. 8 (2022)

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

wheretj = 1 −(ℓj − λs 2 )+
ℓi



2
(ℓ
−
λ
s
)

j
+
 .
;e
tj = − 1l{ℓj > λs 2 } +
2
ℓj
ℓ
1 
j

29 / 31
MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

Concluding remarks
1

Derivation of an adaptive threshholding rule:For λ ≥ 0, solve a regularized minimizaton problem(random{}one) b1−
1/22e ∈ arg min ∈ Ξ ∥X SΞ− Ξ∥F + λ∥Ξ∥1 .
2Ξ∈Mat(m, p ; R)
(∑
)
b
1
/
2
⊤
We have Ξb =Ξ
eS=ℓ (ℓ − λ) u v S 1/2 .
jλ{wherer (ℓj , uj , v j

X S −1/2 .
Obtain SURE

j =1, 2, ..., m

(

)

+ j j

is a system of singular values of[(bλ , Ξ = E SURE ΞbλRΣ ΞSolve the minimization( )problem
bλ =⇒ Ξbb..
b
λ ∈ arg min SURE Ξ
λ
λ≥0

2

}

j

j =1

)]

It is routine to convert this result to case for complex normal
distribution.
30 / 31

Mathematical Optimization and Statistical Theories using Geometric Methods

121

MANOVA model and its canonical mode
Problem set-up
Mean matrix estimation when a covariance is known
Mean matrix estimation when a covarianc matrix is unknown
Concluding remarks

References
1

Candés, E.J., Sing-Long, C.A., and Trzasko, J.D. (2013):
IEEE on Signal Processing 61 4643–4657.

2

Chételat, D. and Wells, M.T. (2012): AOS 40 3137–3160.

3

Efron, B. (2004): JASA 99 619–642.

4

Hansen, N.R. (2018): SPL 135 76–88.

5

Josse, J. and Sardy, S. (2016): Stat. Comput 26 715–724.

6

Li, K., Li, H., Chen, R.H., and Wen, Y.-W. (2022): SIAM JSC
44 A2204-A2225.

7

Mukherjee, A., Chen, K., Wang, N., and Zhu, L. (2015):
Biometrika 102 457–477.

8

Tsukuma, H. and Kubokawa, T. (2015): JMVA 139 312–328.

31 / 31

122

OCAMI Reports Vol. 8 (2022)

Expected Euler characteristic heuristic for smooth Gaussian
random fields with inhomogeneous marginals
Satoshi Kuriki
The Institute of Statistical Mathematics
10-3 Midoricho, Tachikawa, Tokyo 190-8562, Japankuriki@ism.ac.jp

Abstract
Expected Euler characteristic (EC) heuristic is a method for approximating the tailprobability of the maximum of a Gaussian random field. In this talk, we provide an expected Euler characteristic formula for the approximate tail probability and its relative
approximation error when the index set M is a closed manifold and the mean and varianceof the marginal distribution are not necessarily constant. When the variance is constant,[TTA05] proved that the relative approximation error is exponentially small in a general
setting where the index set M is a stratified manifold. When the variance is not constant,
it is shown that only the subset Msupp of M , referred to as the supporting index set,
contributes to the maximum tail probability. The proposed tail probability formula isan integral of the Euler characteristic density over Msupp , and its relative approximation
error is proven to be exponentially small as in the case of constant variance. These results
are generalizations of [KTT22], who addressed a restricted case of finite Karhunen-Loève
expansion by the volume-of-tube method. As an example, the tail probability formulafor the largest eigenvalues of noncentral Wishart matrices Wp (ν, Σ; Φ) and its relative approximation error are obtained. Numerical experience supports the high accuracy of the
expected Euler characteristic formulas regardless of whether the marginals are homogeneous or inhomogeneous.
Keywords: Borel’s inequality, Kac-Rice formula, noncentral Wishart distribution, volumeof-tube method, Weyl’s tube formula.

References[KTT22] Satoshi Kuriki, Akimichi Takemura, and Jonathan E. Taylor, The volume-oftube method for gaussian random fields with inhomogeneous variance, Journal ofMultivariate Analysis 188 (2022), 104819.
[TTA05] Jonathan E. Taylor, Akimichi Takemura, and Robert J. Adler, Validity of the
expected Euler characteristic heuristic, Ann. Probab. 33 (2005), no. 4, 1362–
1396.

Mathematical Optimization and Statistical Theories using Geometric Methods

PATTERN RECOVERY BY SLOPEPIOTR GRACZYK

Abstract
I will present recent results obtained in [1] jointly with M. Bogdan,
X. Dupuis, B. Kolodziejek, T. Skalski, P. Tardivel and M. Wilczyński.
SLOPE is a popular method for dimensionality reduction in the highdimensional regression. Indeed, some regression coefficient estimatesof SLOPE can be null (sparsity) or can be equal in absolute value(clustering). Consequently, SLOPE may eliminate irrelevant predictors
and may identify groups of predictors having the same influence on thevector of responses.
The notion of SLOPE pattern allows to derive theoretical properties
on sparsity and clustering by SLOPE. Specifically, the SLOPE patternof a vector provides: the sign of its components (positive, negative ornull), the clusters (indices of components equal in absolute value) andclusters ranking.
In this research we give a necessary and sufficient condition for
SLOPE pattern recovery of an unknown vector of regression coefficients.

References[1] M. Bogdan, X. Dupuis, P. Graczyk, B. Kolodziejek, T. Skalski, P. Tardivel,
M. Wilczyński, Pattern recovery by SLOPE (2022), arXiv:2203.12086.
Université d’Angers, CNRS, LAREMA, SFR MATHSTIC, F-49000 Angers,
FranceEmail address: graczyk@univ-angers.fr

123

124

OCAMI Reports Vol. 8 (2022)

Pattern Recovery by SLOPE
Piotr Graczyk
LAREMA, Université d’Angers, France

OCAMI Workshop 20 - 21 October 2022
Mathematical optimization and statistical theories using
geometric methods

1/14

Piotr Graczyk

Pattern Recovery by SLOPE

Joint work with:
T.Skalski (Angers and Wroclaw joint PhD)
M. Bogdan, M. Wilczyński (Wroclaw)
B. Kolodziejek (Warsaw)
X. Dupuis, P. Tardivel (Dijon)
U. Schneider (Vienna)

2/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods[1] M. Bogdan, X. Dupuis, P. Graczyk, B. Kolodziejek, T. Skalski,
P. Tardivel, M. Wilczyński, Pattern recovery by SLOPE (2022),
arXiv:2203.12086.
This paper is purely analytical, even if some intuitions and notions aregeometrical.

[2] P. Tardivel, T. Skalski, U. Schneider, P. Graczyk,
The Geometry of Model Recovery by Penalized and Thresholded
Estimators (2022), HAL preprint hal-03262087.
...............................................
A geometrical approach to SLOPE was initiated in[S-T] U. Schneider, P. Tardivel(2020). The Geometry of
Uniqueness, Sparsity and Clustering in Penalized Estimation. arXiv
preprint arXiv:2004.09106, to appear in 2022.

3/14

Piotr Graczyk

Pattern Recovery by SLOPE

Linear regression model

We dispose of n observations of p explicative variables (predictors)
X1 , . . . , Xp and a response variable Y :Yi = β1 xi1 + · · · + βp xip + εi , i = 1, . . . , n.
X = (xij )1≤i≤n,1≤j≤p is the design n × p matrix.
The columns of X correspond to p variables

β = (β1 , . . . , βp ) ∈ Rp unknown regression coefficients.
ε = (ε1 , . . . , εp ) ∈ Rn random noise.

Matrix notation: Y = X β + ε

4/14

Piotr Graczyk

Pattern Recovery by SLOPE

125

126

OCAMI Reports Vol. 8 (2022)

Linear regression Y = X β + ε, X ∈ Rn×p .
Estimator of β?
Classical statistics case: p ≤ n, rankX = p
Ordinary Least Squares estimator:β̂ OLS = arg minb∈Rp ∥Y − Xb∥22

5/14

Piotr Graczyk

Pattern Recovery by SLOPE

Linear regression Y = X β + ε, X ∈ Rn×p .
Estimator of β?
Classical statistics case: p ≤ n, rankX = p
Ordinary Least Squares estimator:β̂ OLS = arg minb∈Rp ∥Y − Xb∥22 = (X ′ X )−1 X ′ Y
Chalenging case: p > nβ̂ OLS is not uniquely determined, so no longer useful
Modern statistics resorts to the penalized least squares estimators:β̂ = arg min ∥Y − Xb∥22 + pen(b),b∈Rpwhere pen is the penalty on the model complexity.

5/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Penalized estimators LASSO and SLOPELASSO (Tibshirani (1996)):

6/14

pen(b) = λ∥b∥1 , λ > 0

Piotr Graczyk

Pattern Recovery by SLOPE

Penalized estimators LASSO and SLOPELASSO (Tibshirani (1996)):pen(b) = λ∥b∥1 , λ > 0

SLOPE (Sorted L One Penalized Estimation)(Bogdan et al. (2015)), defined as
p

β̂

SLOPE

X
1
λi |b|(i) ,= argmin ∥Y − Xb∥22 +b∈Rp 2|i=1 {z }sorted ℓ1 norm

where λ1 > 0, λ1 ≥ . . . ≥ λp ≥ 0 and |b|(1) ≥ . . . ≥ |b|(p) .

6/14

Piotr Graczyk

Pattern Recovery by SLOPE

127

128

OCAMI Reports Vol. 8 (2022)

Penalized estimators LASSO and SLOPELASSO (Tibshirani (1996)):pen(b) = λ∥b∥1 , λ > 0

SLOPE (Sorted L One Penalized Estimation)(Bogdan et al. (2015)), defined as
p

β̂

SLOPE

X
1
2
λi |b|(i) ,= argmin ∥Y − Xb∥2 +b∈Rp 2|i=1 {z }sorted ℓ1 norm

where λ1 > 0, λ1 ≥ . . . ≥ λp ≥ 0 and |b|(1) ≥ . . . ≥ |b|(p) .
When λ1 = . . . = λp > 0 then SLOPE coincides with LASSO.
Our results for SLOPE give a new approach to LASSO.

6/14

Piotr Graczyk

Pattern Recovery by SLOPE

Polyhedral penalties and dimensionality reduction

In case when the penalty function pen is a polyhedral norm

7/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Polyhedral penalties and dimensionality reduction

In case when the penalty function pen is a polyhedral norm
(i.e. the unit ball Bpen (0, 1) ⊂ Rp in the pen norm is a polyhedron)

7/14

Piotr Graczyk

Pattern Recovery by SLOPE

Polyhedral penalties and dimensionality reduction

In case when the penalty function pen is a polyhedral norm
(i.e. the unit ball Bpen (0, 1) ⊂ Rp in the pen norm is a polyhedron)penalized estimators usually possess the dimensionality reduction
properties.

7/14

Piotr Graczyk

Pattern Recovery by SLOPE

129

130

OCAMI Reports Vol. 8 (2022)

Polyhedral penalties and dimensionality reduction

In case when the penalty function pen is a polyhedral norm
(i.e. the unit ball Bpen (0, 1) ⊂ Rp in the pen norm is a polyhedron)penalized estimators usually possess the dimensionality reduction
properties.
It is well known that LASSO estimator has many null components
β̂iLASSO = 0.
Dimensionality reduction property of LASSO consists inelimination of irrelevant predictors Xi .

7/14

Piotr Graczyk

Pattern Recovery by SLOPESLOPE: dimensionality reduction also by clustering
variables
Another important kind of dimensionality reduction consists inclustering (merging, summing) variables with the same values ofregression coefficients:β̂i = β̂j =⇒ Y = ... + β̂i (Xi + Xj ) + ...
LASSO does not have this property!

8/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

SLOPE: dimensionality reduction also by clustering
variables
Another important kind of dimensionality reduction consists inclustering (merging, summing) variables with the same values ofregression coefficients:β̂i = β̂j =⇒ Y = ... + β̂i (Xi + Xj ) + ...
LASSO does not have this property!
Statisticians working with SLOPE observed that many coefficient
regression estimates of SLOPE can be:

8/14

Piotr Graczyk

Pattern Recovery by SLOPESLOPE: dimensionality reduction also by clustering
variables
Another important kind of dimensionality reduction consists inclustering (merging, summing) variables with the same values ofregression coefficients:β̂i = β̂j =⇒ Y = ... + β̂i (Xi + Xj ) + ...
LASSO does not have this property!
Statisticians working with SLOPE observed that many coefficient
regression estimates of SLOPE can be:• equal =⇒ clustering predictors• null =⇒ eliminating irrelevant predictors like LASSO

8/14

Piotr Graczyk

Pattern Recovery by SLOPE

131

132

OCAMI Reports Vol. 8 (2022)

Simulations: n = 100, p = 200, LASSO and SLOPE on R
We simulated Y = X β + ε where ε has iid N(0, 52 ) entries andβ1 = . . . = β30 = 40,
β31 = . . . = β200 = 0.
The rows of the design matrix X are generated as independent binary Markov chains, with
P(Xi1 = 1) = P(Xi1 = −1) = 0.5 and P(Xi(j+1) ̸= Xij ) = 1 − P(Xi(j+1) = Xij ) = 0.0476.
Both LASSO and SLOPE properly estimate at 0 null components of β (not drawn)
Cluster

●
●

40
38

estimator

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

36

42

44

● LASSO● SLOPE

●

36

38

40

42

44

β

9/14

Piotr Graczyk

Pattern Recovery by SLOPE

Main objective of our research

Why/when does SLOPE recover the clusters
and zeros (” SLOPE pattern”) of β?

10/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Main objective of our research

Why/when does SLOPE recover the clusters
and zeros (” SLOPE pattern”) of β?
Explain this phenomenon strictly
mathematically.

10/14

Piotr Graczyk

Pattern Recovery by SLOPE

Main objective of our research

Why/when does SLOPE recover the clusters
and zeros (” SLOPE pattern”) of β?
Explain this phenomenon strictly
mathematically.
Give sufficient and necessary conditions for
SLOPE pattern recovery.

10/14

Piotr Graczyk

Pattern Recovery by SLOPE

133

134

OCAMI Reports Vol. 8 (2022)

Main objective of our research

Why/when does SLOPE recover the clusters
and zeros (” SLOPE pattern”) of β?
Explain this phenomenon strictly
mathematically.
Give sufficient and necessary conditions for
SLOPE pattern recovery.
A by-product: a new and simple mathematical approach to thesequestions for LASSO (huge literature on LASSO is very technical)

10/14

Piotr Graczyk

copy to come from copywriter for review

143

Monographs on Statistics and Applied Probability 143

Statistical Learning with Sparsity

Statistics

Statistical Learningwith Sparsity

Hastie • Tibshirani • Wainwright

K25103

Pattern Recovery by SLOPE

The Lasso andGeneralizations

Trevor Hastie
Robert Tibshirani
Martin Wainwrightw w w. c rc p r e s s . c o m

K25103_cover.indd 1

2/24/15 1:35 PM

11/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

12/14

Piotr Graczyk

Pattern Recovery by SLOPE

Dimensionality reduction by SLOPE
Some coefficient regression estimates of SLOPE can be null or canbe equal in absolute value.

Figure: This figure intuitively illustrates that SLOPE can have some nullcomponents or some components equal in absolute value.
13/14

Piotr Graczyk

Pattern Recovery by SLOPE

135

136

OCAMI Reports Vol. 8 (2022)

Dual penalty norm and dual ballSuppose that pen is a polyhedral norm on Rp .

14/14

Piotr Graczyk

Pattern Recovery by SLOPE

Dual penalty norm and dual ballSuppose that pen is a polyhedral norm on Rp .
Our results show that the dual unit ball B ∗ plays a crucial role instudying penalized estimators rather than B itself.

14/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Dual penalty norm and dual ballSuppose that pen is a polyhedral norm on Rp .
Our results show that the dual unit ball B ∗ plays a crucial role instudying penalized estimators rather than B itself.
Given a norm ∥ · ∥ on Rp , recall that the dual norm ∥ · ∥∗ is definedby∥b∥∗ = max{v ′ b : ∥v ∥ ≤ 1} = ∥b ∗ ∥,i.e. it is the norm of b considered as a linear functional b ∗ .

14/14

Piotr Graczyk

Pattern Recovery by SLOPE

Dual SLOPE norm and dual ballLet Λ = (λ1 , . . . , λp )′ where λ1 > 0 and λ1 ≥ . . . ≥ λp > 0.
The sorted ℓ1 norm is denotedJΛ (b) =pXi=1

λi |b|(i) where |b|(1) ≥ . . . ≥ |b|(p) .

The dual sorted ℓ1 norm is equal to


|b|
|b|
+
|b|
|b|
+
.
.
.
+
|b|
(1)
(1)
(2)
(1)
(p)
,
,...,
.
JΛ∗ (b) = max
λ1
λ1 + λ2
λ1 + . . . + λp
The dual SLOPE ball is defined by
B ∗ = {v ∈ Rp | JΛ∗ (v ) ≤ 1}.
B ∗ is a signed permutahedron in Rp : its vertices are signed
permutations of Λ.
15/14

Piotr Graczyk

Pattern Recovery by SLOPE

137

138

OCAMI Reports Vol. 8 (2022)

p = 3, B ∗ = signed permutahedron

16/14

Piotr Graczyk

Pattern Recovery by SLOPE

Approach of minimization by subdifferential

Let f : Rp → R be a convex function.
The subdifferential ∂f is defined by∂f (b) = {v ∈ Rp : f (z) ≥ f (b) + v ′ (z − b) ∀ z ∈ Rp }Evidently, f attains its minimum at a point b if and only if
0 ∈ ∂f (b)
.

17/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Approach of minimization by subdifferential
Recall the SLOPE minimization problem:
minimize b → f (b) = 12 ∥Y − Xb∥22 + JΛ (b).
It is a particular case of pen–minimization problem
minimize b → f (b) = 12 ∥Y − Xb∥22 + pen(b).
Proposition (Solution of pen-min problem)β̂ is a solution of the pen minimization problem if and only if
X ′ (Y − X β̂) ∈ ∂(pen)(β̂).
Proof. f attains its minimum at a point b if and only if 0 ∈ ∂f (b).
We have∂f (b) = ∂ 12 ∥Y − Xb∥22 + ∂(pen)(b) = {−X ′ (Y − Xb)} + ∂(pen)(b).
The condition 0 ∈ ∂f (b) gives the proposition.
□.
Thus we need to understand ∂(pen).
18/14

Piotr Graczyk

Pattern Recovery by SLOPE

Subdifferential of a norm and the dual ball B ∗Proposition (Subdifferential and the dual ball)(a) The subdifferential of a norm ∥ · ∥ is the following subset of B ∗ :∂∥ · ∥(b) = {v ∈ Rp : ∥v ∥∗ ≤ 1 and v ′ b = ∥b∥}(b) If the norm ∥ · ∥ is polyhedral, then ∂∥ · ∥(b) is a face of B ∗and all faces of B ∗ are subdifferentials of ∥ · ∥.
Proof. (a) is an easy exercice. Both parts are in the book:
HIRIART-URRUTY, J.-B. and LEMARÉCHAL, C. (2004).
Fundamentals of convex analysis. Springer.

19/14

Piotr Graczyk

Pattern Recovery by SLOPE

139

140

OCAMI Reports Vol. 8 (2022)

Set SX ,Λ (Y ) of SLOPE solutions. Uniqueness.
We denote SX ,Λ (Y ) ̸= ∅ the set of SLOPE solutions. It is easy to
see that it is compact. It may be bigger than a singleton.
The unicity has the following geometrical characterization.
Theorem (Uniqueness, [S-T],[2])The solution of the pen-minimization problem is unique for all
Y ∈ Rn if and only if row(X ) does not intersect a face of the dualball B ∗ whose codimension is greater than dim(col(X )).
• Cases in which SX ,Λ (Y ) is not a singleton are very rare.
Indeed, the set of matrices X ∈ Rn×p for which there exists a
Y ∈ Rn where SX ,Λ (Y ) is not a singleton has a null Lebesgue
measure on Rn×p ([S-T])
If ker(X ) = {0}, then SX ,Λ (Y ) consists of one element.

20/14

Piotr Graczyk

Pattern Recovery by SLOPESLOPE pattern and related notions
The SLOPE pattern ( introduced by Schneider and Tardivel(2020)) extracts from a given vector:a) The sign of the components (positive, negative or null),b) The clusters (indices of components equal in absolute value),c) The hierarchy between the clusters.
Definition (SLOPE pattern)
Let b ∈ Rp . The SLOPE pattern of b, patt(b) ∈ Zp , is defined bypatt(b)i = sign(bi ) rank(|b|)i ,i ∈ {1, . . . , p}

where rank(|b|)i ∈ {0, 1, . . . , k}, k is the number of nonzero
distinct values in {|b1 |, . . . , |bp |}.

21/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Exampleb = (4.7, −4.7, 0, 1.8, 4.7, −1.8)′ → patt(b) = (2, −2, 0, 1, 2, −1)′ .
PpSLOPE = patt(Rp ) denotes the set of SLOPE patterns.

22/14

Piotr Graczyk

Pattern Recovery by SLOPE

Identification of patterns as subdifferentials

Theorem (SLOPE pattern= subdifferential(SLOPE pen))Let Λ = (λ1 , . . . , λp )′ where λ1 > . . . > λp > 0 and a, b ∈ Rp . We
have patt(a) = patt(b) if and only if ∂JΛ (a) = ∂JΛ (b).
Proof. A first (involved) proof was given in [S-T]. In [1] we give a
simple proof as a corollary from the (coming below) Proposition on
affine characterization of ∂(JΛ ) for SLOPE.

23/14

Piotr Graczyk

Pattern Recovery by SLOPE

141

142

OCAMI Reports Vol. 8 (2022)

Identification of patterns as subdifferentials
Consequently,
for any polyhedral norm penalty pen, we define in [2]:Definition (Pattern= subdifferential(pen), [2])For a penalized estimator with pen equal to a polyhedral norm, we
say that patt(a) = patt(b) if a and b have the same
subdifferentials: ∂pen(a) = ∂pen(b).

Example. For LASSO, with pen = ∥ · ∥1 , we getpatt(a) = sign(a).
Indeed, the subdifferentials of pen = ∥ · ∥1 (=faces of the unit ballin ∥ · ∥∞ ) are in bijection with the set {−1, 0, 1}p .
24/14

Piotr Graczyk

Pattern Recovery by SLOPE

Pattern recovery
Definition (Pattern recovery)
We say that SLOPE pattern is recovered by the SLOPE estimator
if there exists β̂ ∈ SX ,Λ (Y ) withpatt(β̂) = patt(β).
Example. Let the true β = (5, 5, 2, −5)′ and the SLOPE estimatorβ̂1 = (4, 4, 3, −4)′ .
Then patt(β̂) = patt(β) = (2, 2, 1, −2)′ and we have the pattern
recovery.
If β̂2 = (4.01, 3.99, 3, −4)′ , then patt(β̂) = (4, 2, 1, −3) ̸= patt(β)
and there is no pattern recovery.
However, it is natural to round up (threshold)β̂2 = (4.01, 3.99, 3, −4)′ ≈ (4, 4, 3, −4)′ .
The thresholded estimator β̂2thresh recovers the pattern of β.
25/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Accessibility of a pattern

Not all the patterns can be realized by β̂ when p > n.
Definition (Accessible pattern)
Let X ∈ Rn×p and pen be a polyhedral norm. We say that β ∈ Rp
has an accessible pattern with respect to X and pen, if there existsy ∈ Rn and β̂ ∈ SX ,pen such that patt(β̂) = patt(β).

26/14Piotr Graczyk

Pattern Recovery by SLOPE

Accessibility of a pattern
Proposition (Geometric characterization of accessible patterns, [2])The pattern of β ∈ Rp is accessible with respect to X and pen if
and only ifrow(X ) ∩ ∂(pen)(β) ̸= ∅.
Proof. ( =⇒ ) When the pattern of β is accessible with respect to
X and pen, there exists y ∈ Rn and β̂ ∈ SX ,pen (y ) such that∂(pen)(β̂) = ∂(pen)(β). Because β̂ is a minimizer,
X ′ (y − X β̂) ∈ ∂(pen)(β̂) = ∂(pen)(β), so that, clearly,col(X ′ ) = row(X ) intersects ∂(pen)(β).
( ⇐= ) If row(X ) intersects the face ∂(pen)(β), then there existsz ∈ Rn such that X ′ z ∈ ∂(pen)(β). For y = X β + z, we have
X ′ (y − X β) = X ′ z, so that β ∈ SX ,pen (y ) and patt(β) is
accessible with respect to X and pen.
□
27/14

Piotr Graczyk

Pattern Recovery by SLOPE

143

144

OCAMI Reports Vol. 8 (2022)

n = 2, p = 3: typically, 17 patterns accessible from 147
The Figure is from [S-T]

colourtype

orange
red
black
pink
purplegreen
blueyellowsegments
segments
segments
segments
pointspointspointspointsintersection 6= ∅row(X) ∩ Fw (±(1, 0, 0))row(X) ∩ Fw (±(1, 1, 1))row(X) ∩ Fw (±(0, 0, 1))row(X) ∩ Fw (±(−1, 0, 1))row(X) ∩ Fw (±(2, 0, −1))row(X) ∩ Fw (±(2, 1, 1))row(X) ∩ Fw (±(1, 1, 2))row(X) ∩ Fw (±(−1, 0, 2))

face intersected isometric to

SLOPE models±{5.5} × P(3.5,1.5)
P(5.5,3.5,1.5)±{5.5} × P(3.5,1.5)
P(5.5,3.5) × [−1.5, 1.5]{5.5} × {3.5} × [−1.5, 1.5]{5.5} × P(3.5,1.5){5.5} × P(3.5,1.5){5.5} × {3.5} × [−1.5, 1.5]±(1, 0, 0)±(1, 1, 1)±(0, 0, 1)±(−1, 0, 1)±(2, 0, −1)±(2, 1, 1)±(1, 1, 2)±(−1, 0, 2)

Table 1: Accessible SLOPE models with respect to X =28/14

Piotr Graczyk

8 5
8
10 1.25 −6and w = (5.5, 3.5, 1.5)0 .

Pattern Recovery by SLOPE

Figure 4: Illustration of the sign permutahedron Pw± (in brown) and the plane row(X) (in light blue).
Because rk(X) = 2 and row(X) does not intersect any vertex of Pw± (the faces with codimension equalto 3), the SLOPE estimator β̂w (y) is unique for all values of y ∈ R2 . Colored segments and points
are the intersections between row(X) and the faces of Pw± , determining the accessible SLOPE models
shown in Table 1. For example, m = (2, 1, 1)0 is an accessible SLOPE model, which implies that there
exists y ∈ R2 for which the SLOPE minimizer β̂w (y) satisfies β̂w (y)1 > β̂w (y)2 = β̂w (y)3 > 0. In
addition, since m = (2, 1, 0)0 is not an accessible model, one cannot pick y ∈ R2 for which the SLOPE
minimizer satisfies β̂w (y)1 > β̂w (y)2 > β̂w (y)3 = 0.

SLOPE pattern matrix

In order to characterize the SLOPE pattern recovery, we will needsome more notions related to a pattern M.
14Definition
Let 0 ̸= M = (M1 , . . . , Mp )′ ∈ PpSLOPE with k = ∥M∥∞ .
Pattern matrix: UM ∈ Rp×k is defined as follows(UM )ij = sign(Mi )1(|Mi |=k+1−j) , i ∈ {1, . . . , p}, j ∈ {1, . . . , k}.
Example. Let M = (1, 2, −2, 0, −1)′ . Then |M|↓ = (2, 2, 1, 1, 0)′



0
1
1
0



−1
0
UM = 


0
0
0 −1

29/14

Piotr Graczyk


1
1

U|M|↓ = 
0
0
0


0
0

1

1
0

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

UM Rk+ gives all vectors with pattern M
For k ≥ 1 we denote by Rk+ = {κ ∈ Rk : κ1 > . . . > κk > 0}.
Definition of UM implies that for 0 ̸= M ∈ PpSLOPE andk = ∥M∥∞ , for b ∈ Rp we havepatt(b) = M ⇐⇒ there exists κ ∈ Rk+ such that b = UM κ.
Example. Let M = (1, 2, −2, 0, −1)′ and κ = (κ1 , κ2 )′ . Then




0
1
κ2
1
   κ1 
0

 κ1



−κ1 
−1
0
UM κ = 
=



 κ
0
 0 
0 2
0 −1
−κ2

30/14

Piotr Graczyk

Pattern Recovery by SLOPE

Clustered matrix X̃M and clustered parameter Λ̃MDefinition (Clustered matrix and Λ– parameter)
Let X ∈ Rn×p , Λ = (λ1 , . . . , λp ) where λ1 > · · · > λp > 0.
Clustered matrix: X̃M = XUM .

Clustered parameter: Λ̃M = (U|M|↓ )′ Λ.
Example. Let X = (X1 |X2 |X3 |X4 |X5 ), M = (1, 2, −2, 0, −1)′ andΛ = (λ1 , λ2 , λ3 , λ4 , λ5 )′ where λ1 > λ2 > λ3 > λ4 > λ5 > 0.

λ1 + λ2X̃M = (X2 − X3 |X1 − X5 ) and Λ̃M =.
λ3 + λ4
The clustered design matrix X̃M has only k = 2 columns instead of
p = 5.

31/14

Piotr Graczyk

Pattern Recovery by SLOPE

145

146

OCAMI Reports Vol. 8 (2022)

If patt(β) = M, then X β = XUM κ = X̃M κ for κ ∈ Rk+ . In
particular,(i)null components Mi = 0 lead to discard the column Xi fromthe design matrix X ,(ii)

a cluster K ⊂ {1, . . . , p} of M (component of M equal inabsolute value) leads to replace thePcolumns (Xi )i∈K by one
column equal to the signed sum:sign(Mi )Xi .
i∈K

32/14

Piotr Graczyk

Pattern Recovery by SLOPE

New characterization of ∂(JΛ ) for SLOPE
The next Proposition provides a new and useful formula for thesubdifferential of the sorted ℓ1 norm, via an optimal system ofaffine equations. This representation is crucial for the paper [1].
Proposition (Affine characterization of ∂(JΛ ) for SLOPE)
Let b ∈ Rp and M = patt(b). Then we have the following formula:
n
o
p
∗
′
∂JΛ (b) = v ∈ R : JΛ (v ) ≤ 1 and UM v = Λ̃M .

Moreover,the affine spaceno generated by ∂JΛ (b) equalsp′v ∈ R | UM v = Λ̃M .

′ v = Λ̃
Example. For M = (1, 2, −2, 0, −1)′ the condition UMMmeans
v2 − v3 = λ1 + λ2 ,v1 − v5 = λ3 + λ4 .

This description is much more performant than the hyperplaneequation v ′ M = JΛ (M) that we saw before!

33/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Proof.

Let us prove the inclusion
n
o
p
∗
′
∂JΛ (b) ⊃ v ∈ R : JΛ (v ) ≤ 1 and UM v = Λ̃M′ v = Λ̃ .
Assume that v ∈ Rp satisfies JΛ∗ (v ) ≤ 1 and UMM
To prove that v ∈ ∂JΛ (b) it remains to establish that b ′ v = JΛ (b) .
Since b = UM s, where s ∈ Rk+ , we have′b ′ v = (UM s)′ v = s ′ UMv = s ′ Λ̃M = JΛ (b).

The proof of the other inclusion is also elementary but longer, weomit it.

34/14

Piotr Graczyk

Pattern Recovery by SLOPE

Characterization of pattern recovery by SLOPE

The characterization of pattern recovery by SLOPE given in thenext Theorem is the main mathematical result of article.
The main statistical results of paper [1] are based thoroughly on
this characterization Theorem.

35/14

Piotr Graczyk

Pattern Recovery by SLOPE

147

148

OCAMI Reports Vol. 8 (2022)

Given a SLOPE minimizer β̂ ∈ SX ,γJλ (Y ) for which
patt(β̂) = M ̸= 0, we observe that the following two simple
properties occur:
Dual ball condition: for π = X ′ (y − X β̂), we have JΛ∗ (π) ≤ 1.
(Actually, we know more: π ∈ ∂(JΛ )(M))
Positivity condition: Consider the vector′ X β̂ = X̃ ′ XU s = X̃ ′ X̃ s, where s ∈ Rk+ .
X̃M
M
M
M M
Thus we have ∃s ∈ Rk+′ X β̂ = X̃ ′ X̃ s.
X̃M
M M

Getting rid of β̂ in the two conditions by some simple
algebraic operations, including:• the Moore-Penrose pseudo-inverse A+ of A′ )+ X̃ ′ = X̃ X̃ + , the projector onto the space col(X̃ )• P̃M = (X̃M
M M
MM
we derive the necessity of two conditions of the next Theorem.
It is next easy to show that these two conditions are alsosufficient for the recovery of the pattern M.
36/14

Piotr Graczyk

Pattern Recovery by SLOPE

Theorem (Characterization of SLOPE pattern recovery bypositivity and dual ball conditions)
Let X ∈ Rn×p , 0 ̸= β ∈ Rp , Y = X β + ε for ε ∈ Rn , Λ ∈ Rp+ .
Let M = patt(β) ∈ PpSLOPE and k = ∥M∥∞ . Define′ +π = X ′ (X̃M
) Λ̃M + X ′ (In − P̃M )Y .

There exists β̂ ∈ SX ,Λ (Y ) with patt(β̂) = patt(β) if and only if
the two conditions below hold true:

k+
′
′

there exists s ∈ R such that X̃M Y − Λ̃M = X̃M X̃M s(positivity condition), ∗JΛ (π) ≤ 1(dual ball condition).

If the positivity and ball conditions are satisfied, thenβ̂ = UM s ∈ SX ,Λ (Y ) and π = X ′ (Y − X β̂).
37/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Proof of necessity of two conditions for model recovery.
Let β̂ ∈ SX ,JΛ (Y ) with patt(β̂) = M, i.e. β̂ = UM s, s ∈ Rk+ .
′ X β̂ fromWe have X ′ (Y − X β̂) ∈ ∂JΛ (M). We want to deduce X̃M
this inclusion.
′ , by the affine characterization ofMultiplying it by UM
subdifferential, we get′ (Y − X β̂) = Λ̃′′X̃M
M and X̃M X β̂ = X̃M Y − Λ̃M .
The positivity condition is proven.
′ )+ to the last equality X̃ ′ X β̂ = X̃ ′ Y − Λ̃Apply (X̃M
M and use
M
M
′
+
′
˜
the fact that PM = (X̃M ) X̃M is the projector onto col(X̃M ). We
have X β̂ = X̃M s ∈ col(X̃M ) so that P̃M X β̂ = X β̂. We get
′ + ′
′ +
′ +
(X̃M
) X̃M X β̂ = P̃M Y − (X̃M
) Λ̃M ⇒ X β̂ = P̃M Y − (X̃M) Λ̃M

We insert this formula for X β̂ in′ )+ Λ̃ + X ′ (I − P̃ )Y .
B ∗ ∋ X ′ (Y − X β̂) = X ′ (X̃M
M
M
We proved the dual ball condition.
Piotr Graczyk

38/14

Pattern Recovery by SLOPE

Necessary condition for model recovery: Λ̃M ∈ col(X̃M′ )Observe that the positivity condition:
there exists s ∈ Rk+ such that

implies that the property′ Y − Λ̃ = X̃ ′ X̃ s
X̃M
M
M M

′
Λ̃M ∈ col(X̃M)′ (X̃ ′ )+ Λ̃ = Λ̃ )(or equivalently, the projector X̃M
M
M
M
is necessary for the positivity condition.
′ ) automatically holds when n ≥ k and
The condition Λ̃M ∈ col(X̃M′ ) = Rk .
col(X̃M

39/14

Piotr Graczyk

Pattern Recovery by SLOPE

149

150

OCAMI Reports Vol. 8 (2022)

Essential term X ′ (X̃M′ )+ Λ̃M in the dual ball condition′ )+ Λ̃The first term X ′ (X̃M
M in the expression′ )+ Λ̃ + X ′ (I − P̃ )Yπ = X ′ (X̃MnM
M
is essential for the dual ball condition. Actually, the second term

X ′ (In − P̃M )Y = X ′ (In − P̃M )X β + X ′ (In − P̃M )ε = X ′ (In − P̃M )εwill be shown neglectable, under natural conditions on the (strong)signal β or when n → ∞.

40/14

Piotr Graczyk

Pattern Recovery by SLOPE

Noiseless case

The second term is null in the noiseless case ε = 0.
′ )+ Λ̃ ) ≤ 1
The dual ball condition becomes JΛ∗ (X ′ (X̃M
M

We check that the positivity condition holds for αΛ with Λ′ )verifying the necessary condition Λ̃M ∈ col(X̃M
and α > 0 small enough.

We prove the following characterization of SLOPE pattern recovery
in the noiseless case.

41/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

SLOPE IR ⇐⇒ noiseless pattern recovery
Corollary (SLOPE IR ⇐⇒ pattern recovery for ε = 0)

Consider the noiseless case when ε = 0.
There exists α > 0 such that SLOPE with tuning parameter αΛrecovers the pattern patt(β) = M if and only if′ +′JΛ∗ (X ′ (X̃M
) Λ̃M ) ≤ 1 and Λ̃M ∈ col(X̃M).

Then ∃α0 such that for all 0 < α < α0 , SLOPE with tuning
parameter αΛ recovers the pattern of β.
By analogy to LASSO terminology (Zou, Wainwright, de Geer) we
say that the SLOPE Irrepresentability(IR) Condition holds if′ )+ Λ̃ ) ≤ 1 and Λ̃ ∈ col(X̃ ′ )JΛ∗ (X ′ (X̃M
M
M
M′ )+ Λ̃ ∈ ∂J (M))).
( or equivalently X ′ (X̃M
MΛWhen ker(X̃M ) = {0} then the SLOPE IR condition reads:′JΛ∗ (X ′ X̃M (X̃M
X̃M )−1 Λ̃M ) ≤ 1.

42/14

Piotr Graczyk

Pattern Recovery by SLOPE

Example, p = 2, n ≥ 2
Let X = (X1 |X2 ) ∈ Rn×2 such that


1
0.6
X ′X =.
0.6 1
Let Λ = (4, 2)′ , β = (5, 3)′ , M = patt(β) = (2, 1)′ .
X̃M = X and Λ̃M = Λ.
ker(X̃M ) = {0} and′ +JΛ∗ (X ′ (X̃M
) Λ̃M ) = JΛ∗ (X ′ X (X ′ X )−1 Λ) = JΛ∗ (Λ) = 1 ≤ 1.

The SLOPE irrepresentability condition holds true, so the noiselesspattern recovery holds for for αΛ.
Using R, we see that 0 < α < 0.4 garantees the pattern recovery.

43/14

Piotr Graczyk

Pattern Recovery by SLOPE

151

152

OCAMI Reports Vol. 8 (2022)

Noiseless pattern recovery holds for β = (5, 3)′ , pattern= (2, 1)′

44/14

Piotr Graczyk

Pattern Recovery by SLOPE

IR does not hold for β = (5, 0)′ , pattern = (1, 0)′J(4,2) (X ′ X̃M′+ Λ̃M ) = 6.4/4 > 1

45/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

153

Geometrical meaning of π1 := X ′ (X̃M′ )+ Λ̃M′ )+ Λ̃
Proposition (π1 := X ′ (X̃M
M is a meeting point )′ ).
Suppose that Λ̃M ∈ col(X̃M
Then {π1 } = aff(∂JΛ (M)) ∩ col(X ′ X̃M ).

Proof. We use the Proposition on Affine characterization to π1 .
′ (X̃ ′ )+ is the projection on col(X̃ ′ ) we haveSince X̃M
M
M
′
′
′ +
UM
π1 = X̃M(X̃M) Λ̃M = Λ̃M .

Thus π1 ∈ aff(∂JΛ (M)).
′ )+ ) = col(X̃ ),Moreover, since col((X̃MM′′+
we deduce that π1 ∈ col(X (X̃M ) ) = col(X ′ X̃M ).
We omit the (short) proof of unicity of the meeting point.

46/14

Piotr Graczyk□

Pattern Recovery by SLOPE∗
Back to the Example: β = (5, 0)′ , pattern = (1, 0)′ J(4,2)(π1 ) > 1,
the meeting point π1 is not in the pattern face ∂JΛ (M)

β2
col(X ′ X̃M )π1∂JΛ (M)
β1

47/14

Piotr Graczyk

Pattern Recovery by SLOPE

154

OCAMI Reports Vol. 8 (2022)

Meeting point IR for any polyhedral pen
For SLOPE, the space col(X ′ X̃M ) = X ′ X colUM = X ′ X linCM
where CM = UM Rk+ is the ”pattern set” of all x ∈ Rp with thesame pattern as M, i.e.
∂JΛ (x) = ∂JΛ (M)The ”pattern set” can be defined for any penalty pen.
The meeting point π1 of aff∂pen(x) and XX ′ linCM is well definedfor any penalty pen.
In [2] we conjecture that the condition π1 ∈ ∂pen(x) is equivalent
to the Noiseless pattern recovery for any polyhedral pen.
(proof at finish)

48/14

Piotr Graczyk

Pattern Recovery by SLOPELASSO analogues of our SLOPE characterization Theorem
and our SLOPE IR condition
Consider the LASSO sign recovery (i.e. existence of estimatorβ̂ LASSO such that sign(β̂ LASSO ) = sign(β) = S ∈ {−1, 0, 1}p )The LASSO analogue of our characterization Theorem withpositivity and dual ball conditions is new. In conclusion we get
Corollary (New LASSO Irrepresentability condition )Consider the noiseless case when ε = 0.
There exists λ > 0 such that LASSO with tuning parameter λrecovers sign(β) = S if and only if∥X ′ (X̃S′ )+ 1Rk ∥∞ ≤ 1

and 1Rk ∈ col(X̃S′ ).

Here X̃S′ is the design matrix X signed and reduced according to S.
Example. If S = (1, 0, −1, 0)′ and X = (X1 |X2 |X3 |X4 ), then
X̃S′ = (X1 , −X3 ).
49/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

New and old LASSO IR condition
The two conditions∥X ′ (X̃S′ )+ 1Rk ∥∞ ≤ 1 and 1Rk ∈ col(X̃S′ )equivalent to noiseless LASSO sign recovery are new.
When ker(X̃S ) = {0} then 1k ∈ col(X̃S′ ) occurs and∥X ′ (X̃S′ )+ 1k ∥∞ ≤ 1 is equivalent to∥XI′ XI (XI′ XI )−1 SI ∥∞ ≤ 1
where I = supp(S), I = {1, . . . , p} \ I( MI denotes the submatrix of M obtained by keeping columns
corresponding to indices in I )This latter expression is known in literature as the LASSO
irrepresentability condition (Fuchs, Zhao, Zou, Wainwright, deGeer).
50/14

Piotr Graczyk

Pattern Recovery by SLOPE

Symmetric error. Necessity of the SLOPE IR Condition

Corollary
Let Y = X β + ε where ε and −ε have the same distribution.
′ )+ Λ̃ ) > 1 or Λ ∈′If JΛ∗ (X ′ (X̃MM
M / col(X̃M ) then the probability ofpattern recovery by SLOPE is smaller than 1/2.
For LASSO, a similar result when ker X̃S = {0}, was obtained by
Wainwright (2009).

51/14

Piotr Graczyk

Pattern Recovery by SLOPE

155

156

OCAMI Reports Vol. 8 (2022)

Asymptotic Pattern Recovery (Pattern Consistency) whenε ̸= 0. Open IR Condition.
In order to give a sufficient condition for pattern recovery, we muststrengthen SLOPE IR condition to an Open SLOPE IR condition(this also happens with LASSO)
Recall that our SLOPE IR condition is equivalent to′ +X ′ (X̃M) Λ̃M ∈ ∂JΛ (M)The Open SLOPE IR condition is′ +X ′ (X̃M) Λ̃M ∈ ri(∂JΛ (M))

where ri(F ) is the relative interior of F .

Piotr Graczyk

52/14

Pattern Recovery by SLOPE

Open IR Condition is numerically effective′ )+ Λ̃ ∈ ri(∂J (M))
The Open IR Condition X ′ (X̃MM
Λis equivalent to the following computationally verifiable conditions:( ∗ ′ ′ +′JnΛ (X (X̃M ) Λ̃M ) ≤ 1 and Λ̃M ∈ col(X̃M ),
o
Pi
Pi
′
′
+
i ∈ {1, . . . , p} := ∥M∥∞ .
j=1 |X (X̃M ) Λ̃M |(j) =j=1 λjWe count the number of equalities in p inequalities equivalent toJΛ∗ (b) ≤ 1. Recall thatJΛ∗ (b) = max

53/14



|b|(1) |b|(1) + |b|(2)|b|(1) + . . . + |b|(p)
,
,...,
λ1
λ1 + λ2
λ1 + . . . + λp

Piotr Graczyk

Pattern Recovery by SLOPE.

Mathematical Optimization and Statistical Theories using Geometric Methods

157

Asymptotic Pattern Recovery (Pattern Consistency) whenε ̸= 0: Open IR, big tuning and strong signal are sufficient
SX ,αΛ (Y ) = argminb∈Rp

1∥Y − Xb∥22 + αJΛ (b).
2Theorem (Pattern consistency with X fixed)
Let X ∈ Rn×p , 0 ̸= M ∈ PpSLOPE , and Λ = (λ1 , . . . , λp )′ where
λ1 > . . . > λp > 0. (β (r ) )r ≥1 sequence with pattern M:
β (r ) = UM s (r )

(r )

(r )

with s1 > . . . > sk > 0 and k = ∥M∥∞ ,


r →∞
(r )
(r )
∆r = min1≤i<k si − si+1 −→ ∞. STRONG SIGNALLet Y (r ) = X β (r ) + ε, where ε is an arbitrary vector in Rn . If
αr → ∞, αr /∆r → 0 as r → ∞ and′ )+ Λ̃ ∈ ri(∂J (M)),
X ′ (X̃M
OPEN IRMΛthen ∃ r0 > 0 ∀ r ≥ r0 ∃ β̂ ∈ SX ,αr Λ (Y (r ) ) such that patt(β̂) = M.
54/14

Piotr Graczyk

Pattern Recovery by SLOPE

Pattern consistency with p fixed and n → ∞We suppose:
X = Xn random, satisfying a natural Lindeberg-Feller condition;
an incremental error εn = (ϵ1 , . . . , ϵn )′ , where (ϵi )i are i.i.d.
centered with finite variance;(Xn )n and (ϵn )n are independent.
Theorem (Pattern consistency with n → ∞)Let X ∈ Rn×p such that n1 X ′ X → C almost surely when n → ∞,0 ̸= β ∈ Rp and M = patt(β). If limn→∞ αnn = 0,αnlimn→∞ √= ∞ andn′
CUM (UM
CUM )−1 Λ̃M ∈ ri(∂JΛ (M))

then

P

patt(β̂nSLOPE ) −→ patt(β).
55/14Piotr Graczyk

Pattern Recovery by SLOPE

158

OCAMI Reports Vol. 8 (2022)

Strong Pattern consistency with n → ∞, for all ω
Assume additionnally that the rows of Xn are independent and thateach row of Xn has the same law as ξ, where ξ is a random vector
whose components are linearly independent a.s. and thatE[ξi2 ] < ∞ for i = 1, . . . , p.
Theorem (Strong Pattern consistency with n → ∞)Let X ∈ Rn×p such that n1 X ′ X → C almost surely when n → ∞,0 ̸= β ∈ Rp and M = patt(β). If limn→∞ αnn = 0,
limn→∞ √ αn
= ∞ andn log log(n)′CUM (UM
CUM )−1 Λ̃M ∈ ri(∂JΛ (M))

then

∀ω

patt(β̂nSLOPE ) −→ patt(β).
56/14Piotr Graczyk

Pattern Recovery by SLOPE

Simulation study: example already seen
Consider Y = X β + ε where ε has iid N(0, 52 ) entries andβ1 = . . . = β30 = 40 and β31 = . . . = β200 = 0.
′ X̃ )−1 Λ̃ ∈ ri(∂J (M)).
X ′ X̃M (X̃M
M
M
Λ
′
−1
′
∥XI XI (XI XI ) sign(βI )∥∞ ≤ 1.
Cluster

●
●

40
38

estimator

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

36

42

44

● LASSO● SLOPE

●

36

38

40

42

44

β

57/14

Piotr Graczyk

Pattern Recovery by SLOPE

Mathematical Optimization and Statistical Theories using Geometric Methods

Consider Y = X β + ε where ε has iid N(0, 52 ) entries andβ1 = . . . = β100 = 40 and β101 = . . . = β200 = 0.
′ X̃ )−1 Λ̃ ) > 1.
JΛ∗ (X ′ X̃M (X̃M
M
M
′
′
−1
∥XI XI (XI XI ) sign(βI )∥∞ > 1.
Cluster

100

120

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

● LASSO

60
0

20

40

estimator

80● SLOPE

●
●
●
●
●
●
0

10

20

30

40

β

Theorem [2]. Under the accessibility condition THRESHOLDEDSLOPE asymptotically recovers the SLOPE pattern of β
58/14

Piotr Graczyk

Pattern Recovery by SLOPERESEARCH PROGRAMRESEARCH PROGRAM (planned with H. Ishi, B.
Kolodziejek, H. Nakashima)
Study of Pattern recovery for Graphical SLOPE on
Graphical Gaussian Models
Pattern = clusters of equal terms and blocks of 0’s⇐⇒ Colored Graphical Models

59/14

Piotr Graczyk

Pattern Recovery by SLOPE

159

160

OCAMI Reports Vol. 8 (2022)

Mathematical optimization and statistical theories
using geometric methods
Date: October 20–21, 2022 (Japan Standard Time)
Venue : Academic Extension Center (Osaka Metropolitan University)
Contents : Workshop (Hybrid: physical/virtual)• This workshop is held as a part of OCAMI Joint Usage/Research (JPMXP0619217849)“MEXT Joint Usage/Research Center on Mathematics and Theoretical Physics”• This workshop is also supported by Japan Science and Technology Agency,
CREST“Innovation of Deep Structured Models with Representation of Mathematical Intelligence” in ”Creating information utilization platform by integrating mathematical and information sciences, and development to society”

Organizers: Hideto Nakashima (ISM: hideto (at) ism.ac.jp), Yoshihiko Konno(OMU), Hideyuki Ishi (OMU), Kenji Fukumizu (ISM)
Program• October 20 (Thursday)
13:00–13:50

Shoji Toyota (SOKENDAI)
Invariance Learning based on Label Hierarchy

14:00–14:50

Sho Sonoda (RIKEN AIP)
Ridgelet Transforms for Neural Networks on Manifolds andHilbert Spaces

15:00–15:50

Tomonari Sei (The University of Tokyo)
Ushio Tanaka (Osaka Metropolitan University)
Stein-type distributions on Riemannian manifolds

16:10–17:00

Tomasz Skalski (Wroclaw University of Science and Technology:
LAREMA, University of Angers)
On LASSO and SLOPE estimators and their pattern recovery

17:10–18:00

Carlos Améndola (Technical University of Berlin)
Likelihood geometry of correlation models

Mathematical Optimization and Statistical Theories using Geometric Methods• October 21 (Friday)
9:00– 9:50

Piotr Zwiernik (University of Toronto)
Mixed convex exponential families and locally associated graphicalmodels

11:00–11:50

Koichi Tojo (RIKEN Center for Advanced Intelligence Project)Classiﬁcation problem of invariant q-exponential families on homogeneous spaces

13:50–14:40

Yoshihiko Konno (Osaka Metropolitan University)
Adaptive shrinkage of singular values for a low-rank matrix mean
when a covariance matrix is unknown

14:50–15:40

Satoshi Kuriki (The Institute of Statistical Mathematics)
Expected Euler characteristic heuristic for smooth Gaussian random ﬁelds with inhomogeneous marginals

16:00–16:50

Piotr Graczyk (LAREMA, University of Angers)
Pattern recovery by SLOPE

161