On a convergence property of a geometrical algorithm for statistical manifolds
Shotaro Akaho Hideitsu Hino Noboru Murata


The spin-forbidden vacuum-ultraviolet absorption spectrum of ^14N^15N
    W. Ubachs
    December 30, 2023
=====================================================================



In this paper, we examine a geometrical projection algorithm for statistical inference.
The algorithm is based on Pythagorean relation and it is derivative-free as well as representation-free that is useful in nonparametric cases. 
We derive a bound of learning rate to guarantee local convergence.
In special cases of m-mixture and e-mixture estimation problems,
we calculate specific forms of the bound that can be used easily in practice.



§ INTRODUCTION


Information geometry is a framework to analyze statistical inference
and machine learning<cit.>.
Geometrically, statistical inference and many machine learning
algorithms can be regarded as procedures to
find a projection to a model subspace from a given data point.
In this paper, we focus on an algorithm to find the projection.

Since the projection is given by minimizing a divergence,
a common approach to finding the projection is a gradient-based method<cit.>.
However, such an approach is not applicable in some cases.
For instance, several attempts to extend the information geometrical
framework to nonparametric cases<cit.>,
where we need to consider a function space or each data is represented as
a point process. In such a case, it is difficult to compute
the derivative of divergence that is necessary for gradient-based methods, and
in some cases, it is difficult to deal with the coordinate explicitly.

Takano et al.<cit.> proposed a geometrical algorithm
to find the projection for nonparametric e-mixture distribution,
where the model subspace is spanned by several empirical distributions.
The algorithm that is derived based on the generalized Pythagorean theorem
only depends on the values of divergences.
It is derivative-free as well as representation-free, 
and it can be applicable to many machine learning
algorithms that can be regarded as finding a projection,
but its convergence property has not been analyzed yet.
The first contribution of this paper is to extend the algorithm
to more general cases.
The second contribution is to give a condition for the convergence
of the algorithm, which is given as a bound of learning rate. In the case of
the discrete distribution, we obtain specific forms of the bound that can be used
easily in practice.



§ GEOMETRICAL ALGORITHM




 §.§ Projection in a statistical manifold


Here we briefly review the information
geometry in order to explain the proposed
geometrical algorithm based on generalized Pythagorean theorem<cit.>.

Let (S,g,∇,∇) be a statistical manifold,
where S is a smooth manifold with a Riemannian metric g,
dual affine connections ∇ and ∇.
We consider the case that S is (dually) flat, where there exist
a ∇-affine coordinate θ=(θ_1,…,θ_d)
and a ∇-affine
coordinate η=(η_1,…,η_d).
For a flat manifold, there exist potential functions
ψ(θ) and ϕ(η), and the two coordinates θ
and η are transformed each other by
Legendre transform,

    θ_i = ∂ϕ(η)/∂η_i,   η_i = ∂ψ(θ)/∂θ_i,   ψ(θ)+ϕ(η)-θ_iη_i=0.

A typical example of a flat manifold is an exponential family,
where each member of the manifold is a distribution of a random variable x
with parameter ξ=(ξ_1,…,ξ_d),

    p(x;ξ) = exp(ξ_i F_i(x) - b(ξ)),

where F_i(x) is a sufficient statistics and exp(-b(ξ)) is a normalization
factor.
For the exponential family, there are two dual connections, called e-connection
and m-connection (e: exponential, m: mixture).
If we take the e-connection as the ∇-connection, ∇-affine
coordinate θ is equal to ξ called e-coordinate,
and ∇-affine coordinate called m-coordinate is given by

    ζ_i = ∂ b(ξ)/∂ξ_i =
      E_ξ [F_i(x)] = ∫ F_i(x) p(x;ξ)dx,

where the function b(ξ) becomes a potential function
ψ(θ). Note that if we take the m-connection as ∇,
the relation changes in a dual way, i.e.,
ζ becomes θ and ξ becomes η.

Here, for p∈ S,
we denote the corresponding ∇- and ∇-coordinate
by p and p respectively.
Let us consider a submanifold defined by linear combination of
K points p_1,…, p_K∈ S,

    M = {p |p = w_k p_k,  w_k=1},

where 𝐰=(w_1,…,w_K) is a weight vector whose sum is 1. 
The submanifold M
 is an affine subspace and hence it is called an
∇-autoparallel (or ∇-flat) submanifold.
In particular, if K=2, M is a straight line of ∇-coordinate
that is called ∇-geodesic.

We can also consider another submanifold in the dual coordinate,

    M = {p |p =  w_k p_k,  w_k=1},

which is called a ∇-autoparallel (or ∇-flat) submanifold.
The ∇-geodesic is defined by a straight line of ∇-coordinate.

Now let us define a ∇-projection and a ∇-projection
from a point q∈ S onto
a submanifold M.
The ∇-projection is a point ∈ M such that
∇-geodesic between q and
 is orthogonal to M at  with respect to the Riemannian metric
g_ij(). In the statistical manifold, g_ij is taken as

    g_ij(θ) 
      = ∂^2 ψ(θ)/∂θ_i∂θ_j,

which is equal to Fisher information for exponential family

    g_ij(ξ)=E_ξ[∂log p(x;ξ)/∂ξ_ilog p(x;ξ)/∂ξ_j].

In a similar way, ∇-projection onto a submanifold M
is defined as a point 
so that the ∇-geodesic connecting q and  is
orthogonal to M.


 Let M be a ∇-autoparallel submanifold of a statistical manifold S,
 and the ∇-projection be ∈M from a point q∈ S, then
 for any point p∈M, the following relation holds
 
    D(p,q) = D(,q) + D(p,),

 where D is the canonical divergence defined by
 
    D(p,q)    = ψ(q)+ϕ(p)-pq
       =qp
       =qp.


By exchanging ∇ and ∇, we have a dual relation, i.e,
for a ∇-autoparallel submanifold M,
 the ∇-projection ∈ M from a point q∈ S satisfies
 the relation
 
    D(p,q) =D(,q) +D(p,),

 where p∈ M and
 D is a dual divergence defined by D(p, q) = D(q, p).

From this theorem, we see that a ∇-projection (∇-projection)
onto a ∇-autoparallel (∇-autoparallel respectively)
submanifold is unique
and can be found by minimizing corresponding divergence, i.e.,
the ∇-projection is given by

    = min_p∈M D(p, q)

and the ∇-projection is given by

    = min_p∈ MD(p, q).


For the exponential family (<ref>), taking the e-connection as
∇-connection, the divergence is equal to the Kullback-Leibler divergence,

    D(p, q) = ∫ p(x;ξ(p)) logp(x;ξ(p))/p(x;ξ(q))dx.

If we take the e-connection as ∇ or ∇ connection, the
corresponding projection and autoparallel submanifold is called
an e-projection and an e-autoparallel submanifold,
and similarly, an m-projection and an m-autoparallel submanifold are
defined for the m-connection.



 §.§ Geometrical algorithm for projection


Now we propose a geometrical algorithm to find a ∇-projection
(or ∇-projection) onto a ∇-autoparallel
(and ∇-autoparallel respectively)
submanifold.
To avoid redundant description, we only formulate the ∇-projection onto
a ∇-autoparallel submanifold,
since the dual case can be obtained by only exchanging ∇ and
∇.

In this paper, we impose a restriction on the projection.


 The projection belongs to the convex hull
	    of p_1,…,p_K in
(<ref>) and (<ref>), i.e., all w_k>0.

Although the projection from a point q∈ S does not necessarily
belong to the convex hull of basis vectors in general,
some application such as mixture models that will be
explained
in Sec. <ref> requires this assumption.
We will discuss this restriction in Sec. <ref>.






Suppose a point q∈ S and a ∇-autoparallel
submanifold M⊆ S are given,
let ∈M be a current estimate of the projection ∈M (Fig. <ref>) and
let us define the quantity γ_k,

    γ_k =   D(, q) + D(p_k, ) - D(p_k, q).

From Eq. (<ref>), γ_k=0 if and only if =.
If γ_k < 0, that means  is closer to p_k than ,
w_k should be decreased.
On the other hand, if γ_k > 0, 
 is farther from p_k than , w_k should be increased.

From the consideration above, we can construct the Algorithm <ref>
to find the ∇-projection by optimizing weights {w_k}_
k=1,…,K so that  satisfies the Pythagorean relation (<ref>).



In the algorithm, the function f(γ) is a positive and strictly
       monotonically increasing function
       s.t. f(0)=1, which is
       introduced in order to stabilize the algorithm and a typical
choice of f is a sigmoidal function,

    f(γ) = 2/1+exp(-βγ),   β>0.

A parameter β controls the learning speed and it is related
to convergence characteristics of the algorithm.
AK in the case that m-connection is taken as ∇-connection
was firstly introduced by Takano et al.<cit.> in order to
estimate a nonparametric e-mixture distribution.
The main contribution of this paper is to clarify the relation between
the function f and the convergence property.
In later sections, we prove A2 (and also AK)
is locally stable if the derivative of f at the origin is less than
a certain bound. 
For later theoretical analysis, we show the following Lemma here.

The value γ_k in AK is given by 

    γ_k =  (-)(-p_k),

which means that
γ_k only depends on the points on M, if the true projection  is known.



For any p, r∈M and q∈ S, let us define

    γ(p,q,r)    = D(r,q) + D(p, r) - D(p, q) 
       = qr
         rp
          - {qp}
       =  (r-q)(r-p).

The value γ_k is given by

    γ_k =   γ(p_k, q, ) = 
       (-q)(-p_k)
     
       =  (-
      + -q)(-p_k)
     
       =  (-)(-p_k)
     
          + 
     (-q)
     (-+-p_k)
     
       = (-)(-p_k) 
    - γ(,q,) + γ(p_k, q, ).

From the Pythagorean theorem,
 
    γ(, q, ) = γ(p_k, q, ) = 0,

thus γ_k becomes (<ref>).




§ STABILITY ANALYSIS IN THE CASE OF K=2


We start the analysis from the simplest case of K=2.
As shown later, the case of general K is reduced to this case.
From (<ref>), γ_k is only depends on the points
on the M, and if K=2, M is just a one-dimensional
straight line of η.



 §.§ Behavior of γ_k


In order to derive the condition for convergence, we examine the behavior
of γ_k for a small perturbation.

The weight value w_1 can be regarded as an ∇-coordinate of M,
and let u_1 be the ∇-coordinate that is dual to w_1.
Let 1 be the value of w_1 at the projection point ,
and the current estimation k(=w_k^(t))
is perturbed slightly from k,

    1 = 1+ϵ,   2=2-ϵ=
     (1-1)-ϵ,

then from (<ref>), the value γ_1 is given by

    γ_1 = (w_1()-w_1(p_1))(u_1()-u_1()) =
      (1+ϵ-1)Δ u_1,

where w_1()=1 and w_1(p_1)=1 are the w_1 value at
 and p_1 respectively, and

    Δ u_1 = u_1()-u_1().

When ϵ is small, it can be expanded upto the first order of ϵ,

    Δ u_1 = g(w_1) Δ w_1 + o(ϵ),

where

    Δ w_1 = 1 - 1 = ϵ,

 g(w_1) is Jacobian that is equal to Riemannian metric

    g(w_1) = ∂ u_1/∂ w_1,

and is also obtained by

    g(w_1)
      = E_w_1[(∂log p(x;w_1)/∂ w_1)^2]
      = -E_w_1[∂^2log p(x;w_1)/∂ w_1^2].

As a result, we have

    γ_1 = g(1)(1-1)ϵ + o(ϵ).


Similarly,

    γ_2 = g(1)1ϵ + o(ϵ).




 §.§ The condition for local stability of the A2



In this section, we show the condition for local convergence property
of the A2. Here we call the algorithm is locally stable 
when the amount of sufficiently small perturbation from the optimal
solution is decreased by the algorithm.


A2 is locally stable when it holds

    0 < 2/1(1-1) g(1),

where 1 is the optimal weight.



By the A2, the weight w_1 is updated by

    w_1' = 1f(γ_1),

and its first order expansion is given from Eq. (<ref>) by

    w_1'    = 1(1 + 0γ_1) + o(ϵ) 
       =1 + ϵ + 10 g(1)
     (1-1)ϵ + o(ϵ),

and for w_2,

    w_2'    = 2f(γ_2) 
       = 1-1 - ϵ +
     (1-1)g(1)10ϵ + o(ϵ).


We see that w_1' + w_2'=1+o(ϵ), thus the normalization
procedure is negligible up to the first order of ϵ.

The condition that  is a stable point of the algorithm is
given by 

    |w_1' - 1| < |1-1|=|ϵ|.

From Eq. (<ref>), it is

    |ϵ + 10 g(1)
     (1-1)ϵ| < |ϵ|,

which is equivalent to

    10 g(1)
     (1-1) < 2,

then we have

    0 < 2/1 (1-1) g(1) .



Since the true value  is not known when the algorithm is
applied, we have two approaches.
The one is approximating 1 by the current estimate 1
and use adaptively changing the derivative of f,
which will be examined in sec.<ref>.
The other approach is to use a bound that is independent of 1,
which is available in some special cases.


A2 is locally stable when it holds

    0 < 2/sup_w w (1-w) g(w),

 where we denote w=w_1 for simplicity.




§ SPECIAL CASE: DISCRETE DISTRIBUTION


In the following subsections, 
we give specific forms of the bound df(0)/dγ of Eq. (<ref>)
both for the e-projection and m-projection
by considering a discrete distribution as a specific case.

The discrete distribution is given by

    q(x) =  q_i δ_i(x),    x∈{1,2,…,d},


    q_i = 1,    q_i≥0.

where δ_i(x)=1 when x=i and δ_i(x)=0 otherwise.
We see that the discrete distribution belongs to the exponential family as follows:

    q(x)    = exp( (log q_i) δ_i(x)) 
       = exp( (log q_i) δ_i(x) +
     (log q_d)
     (1-δ_i(x))) 
       = exp(logq_i/q_dδ_i(x)
      + log q_d),

where we have d-1 independent parameters q_1,…,q_d-1 and
one dependent parameter q_d is given by q_d=1- q_i.
 By taking the e-connection as the ∇-connection, q(x) becomes
 the same form as
 Eq. (<ref>) by regarding
 
    F_i(x) = δ_i(x),   ξ_i = logq_i/q_d,
          b(ξ) = -log q_d,    i=1,…,d-1.

 The dual coordinate ζ_i is given by
 
    ζ_i = E_q(x)[F_i(x)] = q_i,    i=1,…,d-1.

The basis vectors in S are denoted by

    p_k(x) =  p_kiδ_i(x),   p_ki = 1,    p_ki≥0.




 §.§ The case of e-projection



First, we take the e-connection as the ∇-connection, then the ∇-projection
onto the ∇-autoparallel submanifold is the e-projection
onto the m-autoparallel submanifold.

The m-autoparallel submanifold spanned by p_k(x) is given by a set
of points whose m-coordinate (<ref>) is given by

    ζ_i =  w_k p_ki,   i=1,…,d-1.

Since ζ_i is the probability value, 
it is equivalent to the mixture distribution of {p_k(x)}

    p(x;𝐰) =  w_k p_k(x),   w_k=1,

where w_k is usually assumed to be positive, which matches
the Assumption <ref>.

The mixture distribution has a lot of applications, in which
complicated distribution is decomposed into sum of
simple component distributions. An important application in the
discrete distribution case is the nonnegative matrix factorization<cit.>,
where a matrix X with nonnegative components is approximated by

    X ≃ D C,

where D and C are also matrices with nonnegative components.
Let Π be the normalization operator by which sum of each column 
components become 1. 
It is known<cit.> that if X=DC,
there exist P and W such that

    Π(X) = P W,

where P and W are matrices with nonnegative components and
sum of each column components is 1.
This means that a set of probability distributions are approximated
by mixture of factor distributions.
In the NMF, D and C are optimized alternatively by fixing the other.
Each optimization problem can be regarded as e-projection to m-autoparallel
manifold.

Note that we consider the e-projection onto an m-autoparallel submanifold
in this paper, since it is natural from the generalized Pythagorean
relation. However, many learning algorithms are
formulated to maximum likelihood
that is equivalent to the m-projection, which is different from
e-projection in the sense that the argument of divergence is reversed.
For the discrete distribution case, 
the m-projection to the m-autoparallel submanifold has a unique solution,
but it does not hold in general.

Now we give a sufficient condition for convergence of the e-projection
onto the m-autoparallel submanifold.

The A2
of the e-projection onto an m-autoparallel submanifold
 for the discrete distribution locally stable
if
 
    0 < 2/∑_i(√(p_1i)-√(p_2i))^2,

 where the right hand side has a constant lower bound √(2).



The m-autoparallel model spanned by K=2 points can be written as
 
    p(x;w) = w p_1(x) + (1-w) p_2(x).

The Riemannian metric at p(x;w) is given by

    g(w)     =
     E_w[(∂log p(x;w)/∂ w)^2]
     
       = ∑_x=1^d 1/p(x;w)(∂
     p(x;w)/∂ w)^2
    
       = ∑_x=1^d  (p_1(x)-p_2(x))^2/p(x;w)
       = ∑_i=1^d  (p_1i-p_2i)^2/w p_1i+(1-w)p_2i.

The denominator of right hand side of Eq. (<ref>) is 

    sup_w w (1-w) 
    ∑_i=1^d  (p_1i-p_2i)^2/w p_1i +
      (1-w) p_2i.

The i-th term

    w (1-w) 
    (p_1i-p_2i)^2/w p_1i +
      (1-w) p_2i

has maximum value (√(p_1i)-√(p_2i))^2
when w=√(p_2i)/(√(p_1i)+√(p_2i)), then
Eq. (<ref>) is bounded from upper by

    ∑_i(√(p_1i)-√(p_2i))^2,

 which is a Hellinger distance between p_1(x) and p_2(x),
and  we obtain the sufficient condition for local stability,

    0 < 2/∑_i(√(p_1i)-√(p_2i))^2,

and the right hand side has a constant lower bound √(2).




 §.§ The case of m-projection


In this subsection, we take the m-connection as the ∇-connection,
then the ∇-projection
onto the ∇-autoparallel submanifold is the m-projection
onto the e-autoparallel submanifold.

The e-autoparallel submanifold spanned by p_k(x) is given by a set
of points whose e-coordinate (<ref>) is given by

    ξ_i =  w_k logp_ki/p_kd
      = ( w_k log p_ki) - log p_kd,   i=1,…,d-1.

Since ξ_i=log (q_i / q_d), it is equivalent to the model
specified by

    p(x;𝐰)∝exp( w_k log p_k(x)),
         w_k=1,

which is a different type of mixture, log linear mixture.

We call this type of mixture as e-mixture, while the mixture
specified by Eq. (<ref>) as m-mixture.
Although the e-mixture has not been studied as intensively as the m-mixture,
it has several good properties such as maximum entropy principle.
Takano et al.<cit.> proposed a nonparametric extension of
the e-mixture and its learning algorithm based on the geometrical algorithm,
which is generalized in this paper.
In the nonparametric e-mixture estimation, the basis distributions are
expressed by the empirical distribution (i.e., sum of delta functions), 
thus the e-mixture of basis distibutions cannot mathematically defined.
Instead,  it is defined by geometrical characteristics
of e-mixture<cit.>. Therefore, it is not possible to obtain
the coordinate explicitly. Because the geometrical algorithm is coordinate-free, and it only requires to calculate divergences, which can be
estimated based on nonparametric entropy estimation<cit.>.
This is a strong motivation to propose the geometrical algorithm.

Here we give a sufficient condition for convergence of the m-projection
onto the e-autoparallel submanifold.

The A2
of the m-projection onto the e-autoparallel submanifold
 for the discrete distribution is locally stable if

    0 ≤32/(max_i logp_1i/p_2i -
      min_i logp_1i/p_2i)^2.



The right hand side does not have a constant lower bound unlike the
e-projection
case, and it is left as an open problem whether there exists any constant
bound.


The e-autoparallel model for K=2 is written as

    p(x;w) = 1/Z(w)exp(w log p_1(x) + (1-w) log p_2(x)
    ),

where w is an e-coordinate, Z(w) is a normalization constant

    Z(w) = ∑_x=1^d exp( w log p_1(x) + (1-w) log p_2(x) ).


Since the discrete distribution

    log p_k(x) = ∑_i=1^d log p_kiδ_i(x),

p(x;w) is written as

    p(x;w)    = 1/Z(w)exp((
     wlog p_1i + (1-w)log p_2i)δ_i(x)) 
       =1/Z(w)exp( (a_i w + b_i)δ_i(x)),

where

    a_i = log (p_1i/p_2i),    b_i = log p_2i,


    Z(w) =  c_i(w),    c_i(w) = exp(a_i w + b_i).

Note that p(i;w) = c_i(w) / Z(w)

The Fisher information for this model can be calculated by

    g(w)    = - E_w[∂^2 log p(x;w)/∂ w^2] 
       = 1/Z(w)∂^2 Z(w)/∂ w^2-
     (1/Z(w)∂ Z(w)/∂ w)^2
     
       = a_i^2 c_i(w)/Z(w) -
     (a_i c_i(w)/Z(w))^2 
       =  a_i^2 p(i;w) - ( a_i p(i;w))^2

The last formula represents the variance of a_i
with respect to the probability weight p(i;w).
From Popoviciu's inequality on variances<cit.>,
g(w) has an upper bound that is independent of w,

    g(w)≤1/4(max_i a_i - min_i a_i)^2.

Since w(1-w)≤ 1/4, we obtain the inequality (<ref>) of the Proposition from Eq. (<ref>).





§ LOCAL STABILITY FOR GENERAL K


We proceed to the general case which include K≥ 2.
First we present the main theorem.



Let w_k^*, k=1,…, K be the optimal parameter.
If the function f satisfies

    0 < 2/K max_k w_k^* (1-w_k^*) g(w_k^*),

 AK is locally stable.


The proof is in the appendix.
Basic strategy of the proof is to show the equivalence between
the AK and a component-wise update algorithm based on
A2. In the process of the proof, a possible refinement of
the AK is also suggested.




































































































§ DISCUSSION


In this section, we will discuss several points related to the proposed framework, (1) relation to gradient descent method, (2) possible refinement
of the algorithm, (3) assumption of the positivity.



 §.§ Relation to gradient descent method


In general optimization problems, a gradient descent method is a simple way to solve the problem. Here, we show that the updates of the 
gradient descent and the proposed algorithm are linearly related. 

The parameter {w_k} should satisfy a constraint 
∑_k=1^K w_k=1. We first replace w_K by 1-∑_k=1^K-1 w_k,
then update w_k for k=1,…, K-1 by

    w_k' = w_k - λ∂ D(q̂, q)/∂ w_k,

and w_K' is obtained by 1-∑_k=1^K-1 w_k'.

The gradient of D(q̂,q) with respect to w_k is given by

    ∂ D(,q)/∂ w_k   =
    ∂/∂ w_k[q]
       =∂/∂ w_k[
    ϕ()-q]
       = ∂η_i(q̂)/∂ w_k∂/∂η_i(q̂)[
    ϕ()-∑_i'=1^d i'qi'].

Since 

    η(q̂) = ∑_k=1^K w_k η(p_k)=∑_k=1^K-1 w_k (η(p_k)-η(p_K)) + η(p_K),
 
and ∂ϕ(η)/∂η_i=θ_i, we have

    ∂ D(q̂,q)/∂ w_k   =
    (-q)(p_k-p_K),

which can be, from Lemma <ref>, represented using γ_k,

    ∂ D(q̂,q)/∂ w_k =γ_K-γ_k.

The amount of update by the gradient descent is

    Δ w_k^G = -λ∂ D(q̂,q)/∂ w_k
    =λ (γ_k - γ_K),
       k=1,…, K-1,

and

    Δ w_K^G = -∑_k=1^K-1Δ w_k^G = 
    λ∑_k=1^K-1(γ_K-γ_k).

On the other hand, the amount of update by the AK is
approximated for small γ_k by

    Δ w_k^A = w_k' - w_k = w_k (f(γ_k)-1) 
    ≃ w_k ∂ f(0)/∂γγ_k.

Since ∑_k w_k γ_k=0, no further normalization is
necessary.
Comparing Eq. (<ref>) and Eq. (<ref>) with Eq. (<ref>),
we see that Δ w_k^G and Δ w_k^A are linearly related.
Unlike the gradient descent method, the proposed framework does not
need explicit calculation of the coodinate.



 §.§ Possible refinement of the algorithm


As explained in Sec. <ref>, the condition for convergence
depends on the true parameter, thus one approach to use the adaptively
change the derivative of f is to replace
the true parameter by its estimate. This approach also requires
to estimate the Fisher information.

Another possibility for the refinement of the algorithm is based
on the analysis in the appendix. It will be shown that the AK
is equivalent to the slower version of the component-wise update
algorithm. More specifically, the amount of the update Δ w_k
is smaller by the factor 1-w_k. Therefore, the update
rule in the AK, w_k'=w_k f(γ_k) can be replaced by

    w_k' = w_k f(γ_k/1-w_k),

which does not change the condition of the convergence.




 §.§ On the assumption of positivity




In Sec. <ref>, we assumed that the projection lies on the
convex hull P⊂M
spanned by the basis vectors. In general, however, the projection
point can be out of P. In such a case, we generalize the problem
to find a point on M that minimizes the divergence,

    = min_p∈ P D(p, q).

When the projection point is out of P, the solution  of this problem
is on the boundary of P and the ∇-geodesic connecting q and 
is not orthogonal to M any more.

The proposed algorithm itself works even in this case, because
the boundary is again a convex hull of a subset of basis vectors.
However, we have to be careful about one thing:
once a certain w_k becomes 0, it cannot take positive value
any longer, which means that if the current estimate reaches to the
boundary that does not include the optimal solution, then the estimator
cannot escape from the boundary.

Without the assumption of w_i>0, the ∇-projection of q to a dual autoparallel submanifold M always
exists uniquely and such a formulation is studied
as e-PCA and m-PCA framework<cit.> or exponential family PCA
in a special case<cit.>.
However, the algorithm proposed in this paper cannot be applied as it is, because it is derived under the assumption. 
One method of update for this general case is as follows:
w_k should be increased for positive γ_k and
should be decreased for negative γ_k. Also, ∑_k w_k=1
should be preserved. Therefore, let 𝒦_P be a set of
indices with positive γ_k and 𝒦_P be
a set of indices with negative γ_k. Then

    w_k'    = w_k + ϵγ_k/∑_l∈𝒦_Pγ_l,
       k∈𝒦_P, 
    
    w_k'    = w_k -  ϵγ_k/∑_l∈𝒦_Nγ_l,
       k∈𝒦_N,

for a learning constant ϵ.
There are several variations of such an update, and we also have to take
care the update does not make w_k out of the domain of w_k.
The investigation of convergence property of the modified algorithm
is left as a future work.



§ CONCLUDING REMARKS

We proposed a geometrical projection algorithm that only requires the
calculation of divergences. We also showed the condition of the local
stability of the algorithm. There are various applications in machine learning
and related areas in which the projection onto an autoparallel submanifold
is needed, and they are left as future works.



§ ACKNOWLEDGEMENT

This work was supported by JSPS KAKENHI Grant Numbers 
17H01793, 19K12111.


















plain









§ APPENDIX  PROOF OF THEOREM <REF>









 §.§ Component-wise algorithm







Suppose the point p_k^† be the point where
the extension line from p_k through  intersects with
the boundary of M (Fig. <ref>).
Let us introduce the notation w_k(p) to specify the point p of M,

    p =  w_k(p) p_k,

for instance w_l(p_k) = δ_k(l) and
w_k()=k. 
Since p_k,  and p_k^†
are on the same line, by taking an appropriate ρ,

    p_k^† = (1-ρ) p_k + ρ

or equivalently,

    w_l(p_k^†) = (1-ρ) w_l(p_k) + ρ w_l().

For the p_k^†, w_k(p_k^†) should
be zero, 

    (1-ρ)w_k(p_k) + ρk = 0,

that is

    ρ = 1/1-k,

and the point p_k^† is given by

    p_k^†=
      1/1-k∑_l klp_l


We can consider the component-wise update
for K=2 by using p_k and p_k^†.




More detailed procedures of (a) and (b) are described later. 
In the algorithm,
the number L controls how each component-wise update converges,
which plays an important role for fast convergence as will be demonstrated
in Sec. <ref>.



If the condition of local stability for K=2 is satisfied for all k,
BK,L is locally stable for any L≥1.



Now let us give the procedures in BK,L.
The K=2 algorithm between p_k and p_k^†, the current
solution  should be represented in the form of

    = ω_k p_k + ω_k^†p_k^†,

where ω_k^†=1-ω_k, and then calculate γ_k
and γ_k^† based on the Pythagorean relation, and then
apply the update (<ref>) and (<ref>).
From Eq. (<ref>),

    = -1-ρ/ρp_k + 1/ρp_k^†
      = kp_k + (1-k)p_k^†,

then the weights for p_k and p_k^† are obtained as
ω_k=k and ω_k^†=1-k respectively.
The update of ω_k is written as

    ω_k' = ω_k f(γ_k),   ω_k^†' = ω_k^† f(γ_k^†),

and then normalization is performed as

    ω_k”= ω_k'/ω_k'+ω_k^†',   ω_k^†”= ω_k^†'/ω_k'+ω_k^†' =
     1-ω_k”.

From Eqs. (<ref>) and (<ref>), we have

    = ω_k p_k + ω_k^†/1-k∑_l
      klp_l.

By updating ω_k and ω_k^† to
ω_k” and ω_k^†” respectively,
then the corresponding update of {w_l}_l=1,…,K is given by

    w_k”   = ω_k”, 
    
     w_l”    = ω_k^†”/1-w_kw_l
     =1-w_k”/1-w_kw_l,    l k.

This update requires to calculate p_k^† (and related values),
which increases the computational complexity.
For later discussions, let us rewrite the algorithm when the amount
of update is sufficiently small.
From the discussion on the analysis of K=2 (Eq.(<ref>)), if the update
Δ w_k = Δω_k = ω_k f(γ_k) - ω_k = w_k f(γ_k) - w_k is sufficiently small,

    Δω_k^† = ω_k^† f(γ_k^†) - ω_k^†≃
       - Δ w_k

holds, where ≃ represents the neglecting higher order terms of
Δ w_k.
By this approximation, the update is simplified as follows:



  
Update rule 1:

    k”   = k + Δ w_k, 
    l”   = 1-k”/1-kl
     = l-Δ w_k/1-kl,    l k.


Note that calculating p_k^† is not
necessary any longer.
Based on the Update rule 1, the algorithm is simplified.



BaK,L behaves similarly to BK,L locally and
it requires smaller computation cost.

If the condition of local stability for K=2 is satisfied for all k,
BaK,L is locally stable for any L≥1.




 §.§ One-side component-wise update


The component-wise update without any approximation
requires to find p_k^†, which
may cause additional complexity compared to the AK.
Here we consider a simpler algorithm:
only the k-th weight is updated with fixing other weights
and normalize all weights, that is,



  
Update rule 2:

    k'    = k + Δ w_k, 
    k”   = k'/k' + ∑_l kl =
     k'/1+Δ w_k, 
    l”   = l/1 + Δ w_k,   l k.


This update does not require the computation of p_k^†.
We examine the relation between Update 1 and 2.
Assuming Δ w_k is sufficiently small, the Update 2 is approximated by

    k”≃k + Δ w_k/1 + Δ w_k≃
      (k + Δ w_k) (1 - Δ w_k) ≃k + (1-k)Δ w_k,

 
    w_l”≃ (1-Δ w_k)l = l - Δ w_k l,    l k,

which means that the Update rule 2 is equivalent to the Update rule 1 where
the learning constant is shortened by a factor 1-k.

Therefore, we see that if the Update rule 1 is locally stable,
the Update rule 2 is also locally stable.

In a similar way with component-wise algorithm, 
we can obtain one-side component-wise algorithm
for general K based on Update rule 2.



Note that updating w_k affects the value of other w_l (l k)
because of the normalization.


If the condition of local stability for K=2 is satisfied for all k,
CK,L is locally stable for any L≥1.





 §.§ Local stability of the AK



Now we are ready to prove the local stability of AK.

The CK,1 is a sequential algorithm, and we can 
construct corresponding “batch” version of
the algorithm.



Since CbK updates the weights by K perturbations,
the condition for local stability is changed.
The following lemma gives a sufficient condition.

CbK is locally stable if f satisfies

    0 < 2/K max_k w_k^* (1-w_k^*) g(w_k^*),

 where w_k^* denotes the optimal value.

This bound is given by multiplying 1/K to (<ref>), but it might
be very strict, because (<ref>) for K=2 is a better bound.
Further, as wee see,
CbK is very similar to CK,1, where only the difference
is whether the former is a simultaneous update and the latter is a sequential update.
It is an open problem to obtain a better bound for K>2.


By the update of w_k, suppose the weight is changed from
 𝐰̂ to 𝐰̂ + Δ𝐰_k,
 where (Δ𝐰_k)_l=0 because of the weight constraint.
If f  satisfies the condition (<ref>), it holds
 
    𝐰̂-𝐰^*^2 > 
      𝐰̂+Δ𝐰_k -𝐰^*^2

for all k.
By multiplying 1/K to the value of
 df(0)/dγ, the change of weights becomes
 𝐰̂ + Δ𝐰_k / K
 in terms of the first order approximation, and
 the simultaneous update of the all weight, it becomes
 𝐰̂ + Δ𝐰_k / K.
Therefore, the new weight satisfies

    𝐰̂+1/KΔ𝐰_k -^2
        =  1/K(𝐰̂+Δ𝐰_k
     -)^2 
       = 1/K^2(𝐰̂+Δ𝐰_k
     -)^2 
       < 1/K^2 K^2 max_k 𝐰̂+Δ𝐰_k
     -^2 
       <   𝐰̂-^2,

 which shows local stability of CbK.


The main theorem <ref> is proved by showing equivalence between
AK and CbK as follows.


When the update amounts of {w_k} are sufficiently small,
The AK is equivalent to CbK.



By the AK, the weights are updated by
 
    w_k' = w_k + Δ w_k,

 
    w_k”   = w_k'/ w_l' = w_k + Δ w_k/
      (w_l + Δ w_l)
       = w_k + Δ w_l/1 + Δ w_l
       ≃ (w_k + Δ w_k) (1 - Δ w_l) 
       ≃ w_k + Δ w_k - Δ w_l w_k.

On the other hand, the update of w_k in the Update rule 2 for small
 Δ w_k is given by (<ref>) and (<ref>).
Therefore, the amount of change of w_l for the update of w_k is given by
 
    Δ w_l(k) = - Δ w_k w_l,

 then summing up them and we have the update of
 of the CbK by
 
    w_k”'    ≃ w_k + (1-w_k)Δ w_k - ∑_l kΔ w_l w_k
      
       = w_k + Δ w_k - Δ w_l w_k,

which coincides the update of (<ref>).




