text:
# Simulation of reproductive risk and emergence of female reproductive cessation

## 1 Introduction

The significance of post-reproductive survival is one of most intriguing evolutionary unanswered questions related to reproduction and aging. A very important example of that is the human menopause. Despite of several theories proposed, the early female reproductive senescence, which has been also observed, although generally less clear-cut, in some other species, for instance, rats, macaques, lions, pilot whales, chimpanzees, elephants  remains one of the great puzzles of evolution theory.

In the specific case of humans, whether this postreproductive life is an adaptive consequence of natural selection, or is a non-adaptive artifact of the rapid increase in longevity over the past few centuries , many theories why the menopause evolved are based on unusual circumstances affecting the human life history (for review, see Ref. and references therein). It has been suggested that menopause could act evolutionarily to protect the ageing woman from the hazards of childbirth. That is, as humans evolved and became able to reach greater ages, there came a point when survival during childbirth began to decline as a function of further ageing and increased frailty . ’Premature’ death of the mother would also put at risk any existing children and their potential for reproduction . At the age when the risk of child-bearing outweighs the benefit of producing progeny, natural selection would favour women who became infertile; thus, the evolutionary pressure for menopause. A different theory is that the human menopause arises from the prolonged infant dependency on the mother, coupled with the risks of late pregnancy and child-bearing, due to the large neonatal brain size. An alternative theory is that menopause enhances fitness by producing post-reproductive grandmothers who can assist their adult daughters. Recently a theoretical model incorporating both of these theories has been sucessful in providing a possible evolutionary explanation of menopause .

Investigations of evolutionary problems by physicists have in fact boomed in the last few years. In particular, they have pioneered the use of techniques derived from the availability of powerful low-cost computers to fulfill the lack of and complement experimentation . Since computer simulations of natural systems can provide much insight into their fundamental mechanisms, they can be used to test theoretical ideas that could be otherwise viewed as too vague to deserve the status of scientific knowledge . Among the many computer models introduced to describe the evolution of populations , the Penna bit-string model , the Redfield model  and Stauffer model  have been stood out for predicting many phenomena in population dynamics.

The sexual version of the Penna model for mutation accumulation theory was first introduced by Bernardes , followed by Stauffer et al.  who adopted a slightly different strategy. We are going to describe the second one. The genome of each (diploid) organism is represented by two computer words. In each word, a bit set to $`1`$ at a position (“locus”) corresponds to a deleterious mutation - a “perfect” strand would be composed solely of zeros. The effect of this mutation may be felt by the individual at all ages equal to or above the numerical order of that locus in the world. As an example, a bit set to $`1`$ at the second position of one of the bit-strings means that a harmful effect may become present in the life history of the organism to which it corresponds after it has lived for two time periods (“years”). In order to count the accumulated number of mutations and compare it with a specified threshold $`T`$, it is necessary to distinguish between recessive and dominant mutations. A mutation is counted if two bits set to $`1`$ at the same position in both strings or if it appears in only one of the bit-strings but at a dominant position. The dominant positions are randomly chosen at the beginning of the simulation and are the same for all the individuals. Reproduction is modeled by the introduction of new genomes in the population. Each female and male becomes able to reproduce after having reached a minimum age, after which a female randomly chooses an able male to breed and it generates a fixed number of offspring at the completion of each period of life. To construct one offspring genome so, each string of the mother genome is cut at a randomly selected position, the same for both strings, and the left part of one is combined with the right part of the other, thus generating two new combinations of the original genes. Finally, $`m_f`$ deleterious mutations are randomly introduced, and then the selection of one of these completes the formation of the haploid gamete from the mother. The same process occurs with the male’s genome, producing the male gamete with $`m_m`$ deleterious mutations. These two resulting bit-string form the offspring genome. The gender of the newborn is then randomly selected, with equal probability for each sex.

The Redfield’s model is an elegant model requiring much less computer time than the Penna model, but having no age structure. It is not a population dynamics model following the lifetime of each individual, but only simulates their probabilities to survive up to reproduction. In Redfield’s program the population is characterized by a distribution of mutations $`P(m)`$, $`m=0,1,\mathrm{}`$, giving the probability that an individual has $`m`$ genetic diseases in the genome. Darwinian selection of the fittest then transforms this $`P(m)`$ into a survivor distribution $`L(m)(1s){}_{}{}^{m}P(m)`$ giving the probability that a survivor has $`m`$ deleterious mutations in the genome, where $`s=0.1`$ is a selection coefficient. In the sexual case $`n`$ new hereditary mutations happen according to the Poisson distribution $`\mu {}_{}{}^{n}\mathrm{exp}(\mu )/n!`$, where $`\mu 1`$ corresponds to the mutation rate per genome per generation. In the sexual case, after selection transformed the progeny $`P(m)`$ distribution into the survivor distribution $`L(m)`$, mutations of a rate $`\mu `$ produce the female distribution $`F(m)`$ and those of a rate $`\alpha \mu `$ give the male distribution $`M(m)`$. Then male gametes (sperm cells) are produced containing half of the male genome; thus, their mutation number $`m_m`$ is roughly half of the number of the number of mutations in the father’s genome. Analogously, the number $`m_f`$ of mutations in the female gametes (egg cells) is roughly half of the mother’s number of mutations. The fusion of two gametes adds these two numbers, $`m=m{}_{m}{}^{}+m_f`$, to produce the mutations in the progeny distribution $`P(m)`$. The cycle is repeated until changes in each distribution become negligible. In fact, the mutations in each gamete are not exactly half of the parent but follow a binomial distribution, which simulates the random selection of the transmitted half of the genome via the process of meiosis, crossover, and mitosis.

Of particular interest here is the Stauffer model, which although being a severe abstraction of biology, is particularly simple and had early successes in reproducing observed features of real populations, such as the Gompertz law , the catastrophic senescence of Pacific Salmon , the vanishing of cod in the northwest Atlantic through over-fishing , and the social needs of humans for a minimum population size .

The present work can represent a step in the direction of understanding more the evolutionary origin of menopause. In this way, the simulations reported here test recent assumptions proposed in a work using the Penna model, in which menopause was found to self-organize from a combination of two effects: a reproductive risk of the mother that increases with the number of active mutations and thus with age; and child care, in the sense that young children die if their mother dies .

## 2 Model

In the sexual version of the Stauffer model , each individual of the population, which consists of males and females, is genetically represented by two integers, the minimum reproduction age $`a_m`$ and the genetic death age $`a_d`$.

Individuals die with certainty if their age reaches $`a_d`$, and they die earlier with probability $`V=N(t)/N_{max}`$ at every iteration, where $`N(t)`$ is the actual population size and $`N_{max}`$ a so-called carrying capacity of the environment taking into account the limitations of food and space (Verhulst factor). For computer simulations this logistic factor has the benefit of limiting the size of the population to be dealt with. At every time step and for each individual a random number between zero and $`1`$ is generated and compared with $`V`$: if this number is smaller than $`V`$ the individual dies, independently of its age and genetic death age $`a_d`$.

If the female succeeds in surviving until the minimum reproduction age $`a_m`$, it generates, with probability $`p_b`$, $`b`$ offspring every iteration until death. The female randomly chooses a male to mate, with age also greater or equal to $`a_m`$. The offspring randomly inherits the $`a_m`$ and $`a_d`$ values from one of the parents, independently, apart from random changes (mutations) by $`\pm 1`$. The sex of the baby is randomly chosen, each one with probability $`50\%`$. This procedure is repeated for each of the $`b`$ offspring.

The probability $`p_b(i)`$ for individual $`i`$ to get $`b`$ children takes into account the well known (see e.g. ) trade-off between fecundity and longevity and is the lower the larger the difference $`a_da_m`$ is:

$$p_b(i)1/(a_d(i)a_m(i))$$

where we require $`0a_m<a_da_{max}`$. We start our simulations with $`a_m=1`$, $`a_d=16`$ (or $`a_d=16`$) and after thousand time steps the model shows an emergence of senescence and a typical age of death .

In order to simulate the evolution of populations with females subjected only to reproductive risk $`R`$, we now introduce the following ingredients in the model:

* Reproductive Risk: At the moment of giving birth, we calculate the reproductive risk of a female to die. This is done through the expression:

$$Risk=\mathrm{min}[1,\frac{\beta (a(i)a_m(i))}{a_d(i)a_m(i)}]$$

where $`a(i)`$ is the actual female age and $`\beta `$ is a predefined proportionality factor, which can reduce or increase the whole risk function. A random number between zero $`r[0,1]`$ is drawn, and then the females dies if and only if $`r<Risk`$.
* Age of cessation of female reproduction $`a_c`$, introduced here as a new variable. In the usual model all individuals can reproduce at every period of life until death. Now, for females, we define a maximum age of reproduction $`a_c`$, which is equal to $`a_d`$ at the beginning of the simulation. It means that - only at the beginning of the simulation \- females can reproduce until the end of their lives (there is no menopause at the beginning). When a female with a given value of $`a_c`$ gives birth to a daughter, the daughter’s value of $`a_c`$ is mutated randomly to $`a_c\pm 1`$.

## 3 Results

To obtain all curves presented here we simulated $`20`$ different populations (samples) during $`10,000,000`$ time steps. The parameters of the simulations are:

* Initial population = $`30,000`$ (half for each sex);
* Maximum population size $`N_{max}=1,500,000`$;
* Number of offspring $`b=3`$;

Simulations with less samples or shorter times sometimes gave misleading stability. In Figure 1, which corresponds to the standard results obtained using the original sexual Stauffer model without reproductive risk, we present the time evolution of the average value of the minimum age of reproduction $`a_m`$, the genetic death age $`a_d`$ and the maximum age of reproduction $`a_c`$ (top to bottom, respectively). From this figure we can observe that in order to maximize the probability $`p_b`$, which happens when the difference $`(a_da_m)`$ is minimized, $`a_d`$ ($`\overline{a_d}6`$) stays a value not much bigger than $`a_m`$ ($`\overline{a_m}2`$). If the death age $`a_d`$ would evolve to a very large value when compared to the minimum age of reproduction $`a_m`$, then the probability to reproduce $`p_b`$ would be very small, because $`p_b`$ decreases for large $`(a_da_m)`$. As well, we also can see in this figure that the menopause age becomes bigger than the genetic death age, since $`a_d`$ can not be greater than $`a_{max}`$ and there isn’t any restriction in the maximum value of $`a_c`$, in such a way that the maximum age of reproduction $`a_c`$ is allowed to reach values bigger than the death age $`a_d`$. This means only that the females can reproduce until the end of their lifes, i.e., the reproductive cessation in females can not be observed in this old version without reproductive risk.

On the other hand, when the reproductive risk is introduced, as can be seen in Figure 2 (with $`a_{max}=1200`$), the menopause age $`a_c`$ stays below the death age $`a_d`$. From this result we can notice that there is a compromise between the lower reproductive risk $`Risk`$ and the probability to reproduce $`p_b`$, which is related to the value of the difference $`(a_da_m)`$ : if $`(a_da_m)`$ reaches a small value, i.e., $`a_d`$ not much bigger than $`a_m`$, a high reproduction rate dominates and a strong selection pressure caused by the high reproductive risk could drive the population to extinction. In this way, in order to soften the effects from the reproductive risk and to warrant the survival of the population, $`a_d`$ evolves to values much bigger than $`a_m`$.

This particular feature that $`a_d`$ plays in the model with or without reproductive risk have been obtained from simulations when two different maximum values $`a_{max}`$ \- the maximum value which $`a_d`$ can reach in the simulation - were used. When we consider $`a_{max}=200`$ (see Fig. 3) the difference $`(a_da_m)`$ is still big enough to produce a strong selection pressure against the females to warrant the survival of the population. In other words, the high reproductive risk drives the females to remain reproducing as long as they live, it means that they will reproduce until the end of their lifes. However, if $`a_d`$ is allowed to reach a very high value (see Fig. 2, where $`a_{max}=1200`$), the opposite occurs and the cessation of reproduction of the females appears, since now the females do not suffer very strong effects from the reproductive risk due the large difference $`(a_da_m)`$. Given that the reproductive risk also increases with the age, it would seem to make sense for the females to cease reproducing at a certain age $`a_c`$ below the death age $`a_d`$ (but big enough to guarantee the survival of the population) in order to avoid completely the effects from the reproductive risk, since when the females stop to reproduce there will be no reproductive risk anymore.

The emergence of female menopause age obtained here (Fig. 2), using only a reproductive risk in giving birth that increases with age, was found with the Penna bit-string model  only when in addition to some kind of reproductive risk also maternal care was considered. Menopause was not observed in  for systems with only maternal care or with only reproductive risk. Furthermore, only $`20\%`$ of the fertile female population had post-reproductive life and it was possible only when a very high period of maternal care was used, which represented almost $`24\%`$ of the maximum age of life in their simulations . Their simulations had been performed under those restrictions suggested $`40`$ years ago by Williams that menopause “may have arisen as a reproductive adaptation to a life-cycle characterized by senescence, unusual hazards in pregnancy and childbirth, and a long period of juvenile dependence” . While the usual justification for the introduction of maternal care on the menopause restricts it to the human being, our simulations also to extend to other species, we don’t take into account maternal care in our simulations, and nevertheless found menopause or analogs. The human being seems to be the only species to which kin assist in care and provisioning of young, since orphans regularly can survive even in adverse environments , our simulations suggest that this special human feature is not necessary for the emergence of menopause or analogs.

Similar results have been also observed when using different functions for the reproductive risk, however in all cases simulated the emergence of menopause and its analogs depends on the suitable choice of parameters within the same model. Many simulations were performed in order to avoid to impose a maximum value for the death age, $`a_{max}`$, but without success: without this condition it was observed that the death age $`a_d`$ evolves to a fixed value below the cessation age of reproduction $`a_c`$, which seems to drift towards infinity.

In addition, specifically concerning the constant $`\beta `$ which can reduce or increase the whole risk function, we didn’t manipulate it adequately in order to reproduce real problems related to the risks of childbearing the case of humans only. Instead, we had focused on the minimum value of $`\beta `$ for which would be possible to notice the emergence of menopause using the simple Stauffer model. Different simulations were made, as it can be seen in Table 1, that show stationary values of $`a_m`$, $`a_d`$ and $`a_c`$ averaged over a prespecified number of time steps.

| $`\beta `$ | $`\overline{a_m(i)}`$ | $`\overline{a_d(i)}`$ | $`a_c(i)`$ |
| --- | --- | --- | --- |
| 0.00 | 1.79 | 6.80 | $`\mathrm{}`$ |
| 0.25 | 1.82 | 7.03 | $`\mathrm{}`$ |
| 0.50 | 1.85 | 7.33 | $`\mathrm{}`$ |
| 0.75 | 1.89 | 7.75 | $`\mathrm{}`$ |
| 1.00 | 1.95 | 8.34 | $`\mathrm{}`$ |
| 1.25 | 2.02 | 9.43 | $`\mathrm{}`$ |
| 1.50 | 2.15 | 11.88 | $`\mathrm{}`$ |
| 1.75 | 8.50 | 1182.81 | 633.69 |
| 2.00 | 11.18 | 1181.78 | 701.96 |
| 1.75 | 4.67 | 191.56 | 235.23 |

*Table I.* $`a_m`$, $`a_d`$ and $`a_c`$ averaged over the last time steps for different values of $`\beta `$ and using $`a_{max}=1200`$. In the last row, $`a_{max}=200`$.

## 4 Conclusions

This paper reports on a simulation of populations under reproductive risk, in the framework of a simple model for biological aging. We show that a sexual population presents cessation of female reproduction as a functional characteristic of the system evolving under only reproductive risk. Thus simple explanations for menopause might be too simple since its existence in these simulations depends on the numerical values of model parameters.

Acknowledgments

I would like to thank Suzana Moss and Dietrich Stauffer for discussions and a critical reading of the manuscript, and DAAD for financial support.


format: text

text:
# Paramagnetic-ferromagnetic transition in a double-exchange model

## I Introduction

The double-exchange (DE) model  is one of the basic ones in the theory of magnetism. Magnetic ordering appears in this model due to exchange coupling between the core spins and the conduction electrons. The Hamiltonian of the model is

$`H={\displaystyle \underset{nn^{}\alpha }{}}t_{nn^{}}c_{n\alpha }^{}c_{n^{}\alpha }J{\displaystyle \underset{n\alpha \beta }{}}𝐒_n\sigma _{\alpha \beta }c_{n\alpha }^{}c_{n\beta },`$ (1)
where $`c`$ and $`c^{}`$ are the electrons annihilation and creation operators, $`𝐒_n`$ is the operator of a core spin, $`t_{nn^{}}`$ is the electron hopping, $`J`$ is the exchange coupling between a core spin and n electrons, $`\widehat{\sigma }`$ is the vector of the Pauli matrices, and $`\alpha ,\beta `$ are spin indices.

The model (the core spins being treated as classical vectors) began to be studied several decades ago . Because of a recent general interest in manganites, it reappeared in the focus of attention (see reviews  and references therein). In particular, application of the Dynamical Mean Field Approximation (DMFA)  substantially advanced our understanding of the properties of the DE model . Most of the papers dealing with the model, starting from classical paper by De Gennes , considered the DE Hamiltonian with infinite exchange (and with the addition of the antiferromagnetic superexchange, which is crucial for the explanation of magnetic properties of manganites). However lately many properties of the model for arbitrary strength of the exchange were analyzed .

In this paper we calculate the temperature of a paramagnetic-ferromagnetic transition $`T_c`$ in a double-exchange model for arbitrary electron dispersion law, concentration and relation between the exchange coupling and the electron band width by formulating and solving the DMFA equations. We treat the core spins as classical vectors.

## II DMFA equations and $`T_c`$

Spins being considered as classical vectors $`𝐒_n=𝐦_n`$ (with the normalization $`|𝐦|^2=1`$), the problem becomes a single electron one. The DE Hamiltonian in a single electron representation can be presented as

$$H_{nn^{}}=t_{nn^{}}J𝐦_n\sigma \delta _{nn^{}}.$$
(2)
Let us introduce Green’s function

$`\widehat{G}(E)=(EH)^1`$ (3)
and local Green’s function

$`\widehat{G}_{\mathrm{loc}}(E)=\widehat{G}_{nn}(E).`$ (4)
In the last Equation, the averaging is with respect to random configurations of the core spins. In the framework of the DMFA approach to the problem (see  and references therein) the local Green’s function is expressed through the the local self-energy $`\widehat{\mathrm{\Sigma }}`$ by the equation

$`\widehat{G}_{\mathrm{loc}}(E)=g_0\left(E\widehat{\mathrm{\Sigma }}(E)\right),`$ (5)
where

$`g_0(E)={\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}\left(Et_𝐤\right)^1`$ (6)
is the bare (in the absence of the exchange interaction) local Green’s function. The self-energy satisfies equation

$`\widehat{G}_{\mathrm{loc}}(E)={\displaystyle \frac{1}{\widehat{G}_{\mathrm{loc}}^1(E)+\widehat{\mathrm{\Sigma }}(E)+J𝐦\widehat{\sigma }}},`$ (7)
where $`X(𝐦)X(𝐦)P(𝐦)`$, and $`P(𝐦)`$ is a probability of a given spin orientation (one-site probability). The quantities $`\widehat{G}`$ and $`\widehat{\mathrm{\Sigma }}`$ are $`2\times 2`$ matrices in spin space.

The difference between Eq. (7) and similar equation appearing in the Coherent Potential Approximation , is due to the fact that in our problem the disorder is annealed (and not quenched). That is the probability of a given (random) core spin configuration is determined by the energy of electron subsystem interacting with this configuration.

To understand the DMFA assumption allowing to find the probability $`P(𝐦)`$ self-consistently with the solution of Eq. (7), notice that Eq. (7) reduces the problem of electron scattering by many spins, to scattering by a single spin with the effective potential $`V=J𝐦\sigma \widehat{\mathrm{\Sigma }}`$. The spin is embedded in an effective medium, described by the Hamiltonian $`t_𝐤+\widehat{\mathrm{\Sigma }}`$, and, hence, by the local Green’s function $`\widehat{G}_{\mathrm{loc}}`$. Consider the change in the number of electron states in such system due to a single spin. We get

$`\mathrm{\Delta }D(E,𝐦)={\displaystyle \frac{1}{\pi }}\mathrm{Im}\mathrm{ln}\mathrm{det}\left[1+\left(J𝐦\widehat{\sigma }+\widehat{\mathrm{\Sigma }}\right)\widehat{G}_{\mathrm{loc}}\right],`$ (8)
where the argument of both $`G_{\mathrm{loc}}`$ and $`\mathrm{\Sigma }`$ is $`E+i0`$. So the change in thermodynamic potential is

$`\mathrm{\Delta }\mathrm{\Omega }(𝐦)={\displaystyle f(E)\mathrm{\Delta }D(E,𝐦)𝑑E},`$ (9)
where $`f(E)`$ is the Fermi function. The DMFA approximation for the one-site probability $`P(𝐦)`$ is:

$`P(𝐦)\mathrm{exp}\left[\beta \mathrm{\Delta }\mathrm{\Omega }(𝐦)\right].`$ (10)
Eqs. (7) and (10) are the system of non-linear (integral) equations. However, the system is grossly simplified in paramagnetic (PM) phase, where the macroscopic magnetization $`𝐌=0`$. In this case $`P(𝐦)=\mathrm{const}`$, the averaging in Eq. (7) can be performed explicitly, and we obtain

$`g(E)={\displaystyle \frac{1}{2}}{\displaystyle \underset{(\pm )}{}}{\displaystyle \frac{1}{g^1(E)+\mathrm{\Sigma }(E)\pm J}},`$ (11)
where $`\widehat{\mathrm{\Sigma }}=\mathrm{\Sigma }\widehat{1}`$, $`\widehat{G}=g\widehat{1}=g_0(E\mathrm{\Sigma })\widehat{1}`$, where $`\widehat{1}`$ is a unity matrix.

In the ferromagnetic (FM) phase near the Curie temperature, Eqs. (7) and (10) can be linearized  with respect to $`𝐌`$. Thus we reduce the DMFA equations to a traditional MF equation

$`P(𝐦)\mathrm{exp}\left(3\beta T_c𝐌𝐦\right).`$ (12)
The parameter $`T_c`$ is formally introduced as a coefficient in the expansion of $`\mathrm{\Delta }\mathrm{\Omega }(𝐦)`$ with respect to $`𝐌`$. Non-trivial solution of the MF equation

$$𝐌=𝐦$$
(13)
can exist only for $`T<T_c`$, hence $`T_c`$ is the ferromagnetic transition temperature. In all cases $`T_c`$ turns out to be much less the chemical potential, so we can consider electron gas as degenerate. In this case

$$_{\mathrm{}}^{\mathrm{}}f(E)\mathrm{}𝑑E_{\mathrm{}}^\mu \mathrm{}𝑑E,$$
(14)
where the chemical potential (in paramagnetic phase) is found from the equation

$$n=\frac{2}{\pi }_{\mathrm{}}^\mu \mathrm{Im}g𝑑E,$$
(15)
and $`n`$ is the number of electrons per site.

For the $`T_c`$, after straightforward, though lengthy algebra, we obtain

$`T_c={\displaystyle \frac{2J^2}{3\pi }}{\displaystyle _{\mathrm{}}^\mu }\text{Im}\left[{\displaystyle \frac{g}{\frac{(\mathrm{\Sigma }g^{}g)(1+\mathrm{\Sigma }g)}{g^{}+g^2}\frac{2J^2g}{3}}}\right]𝑑E,`$ (16)
where $`\mathrm{\Sigma }`$ and $`g`$ are determined by the properties of the system in the PM phase; they are the solutions of Eqs. (5) and (11); $`g^{}\frac{dg_0(\epsilon )}{d\epsilon }|_{\epsilon =E\mathrm{\Sigma }}`$. Eq. (16) is the main result of our paper.

In the case of weak exchange $`J0`$ Eq. (16) takes the form

$`T_c={\displaystyle \frac{2J^2}{3}}\left[N_0(\mu ){\displaystyle \frac{1}{\pi }}{\displaystyle _{\mathrm{}}^\mu }\text{Im}g_0^2(E)\right]dE.`$ (17)
In fact, this equation is the MF approximation  for the RKKY Hamiltonian , to which the original Hamiltonian (1) can be reduced in the case $`JW`$.

In the opposite case $`J\mathrm{}`$ we have to consider only one spin sub-band. (Because of the symmetry of the problem with respect to transformation $`n2n`$, everywhere further on we consider only the case $`n1`$.) After we shift $`E`$, $`\mathrm{\Sigma }`$ and $`\mu `$ by $`J`$, we obtain

$`T_c={\displaystyle \frac{2}{\pi }}{\displaystyle _{\mathrm{}}^\mu }\text{Im}\left[{\displaystyle \frac{g^{}+g^2}{g^{}2g^2}}\right]𝑑E,`$ (18)
where $`g`$ is found from the equation

$`g(E)={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{g^1(E)+\mathrm{\Sigma }(E)}}.`$ (19)
## III $`T_c`$ for semi-circular DOS

Let us apply Eq. (16) to the case of semi-circular (SC) bare density of states (DOS)

$$N_0(\epsilon )=\frac{2}{W}\sqrt{\left(\frac{E}{W}\right)^21}.$$
(20)
Then

$`g_0(E)={\displaystyle \frac{N_0(\epsilon )d\epsilon }{E\epsilon }}={\displaystyle \frac{2}{W}}\left[{\displaystyle \frac{E}{W}}\sqrt{\left({\displaystyle \frac{E}{W}}\right)^21}\right].`$ (21)
From Eq. (5) follows

$$\mathrm{\Sigma }=EW^2g/4g^1,$$
(22)
and Eq. (11) can be written as

$`g={\displaystyle \frac{1}{2}}{\displaystyle \underset{(\pm )}{}}{\displaystyle \frac{1}{EW^2g/4\pm J}}.`$ (23)
Eq. (16) takes the form

$`T_c={\displaystyle \frac{J^2W^2}{6\pi }}{\displaystyle _{\mathrm{}}^\mu }`$ (24)
$`\text{Im}\left[{\displaystyle \frac{g^2}{\left(E\frac{W^2g}{4}\right)\left(E\frac{W^2g}{2}\right)\frac{J^2W^2}{6}g^2}}\right]dE.`$ (25)
In the case of weak exchange $`J/Wn^{1/3}`$, we obtain

$`T_c={\displaystyle \frac{4J^2}{3\pi ^2W}}(4y^21)\sqrt{1y^2},`$ (26)
where

$$n=n_0(y)1+\frac{2}{\pi }\left(\mathrm{sin}^1y+y\sqrt{1y^2}\right).$$
(27)
In the case of strong exchange $`J/Wn^{1/3}`$, Eq. (18) takes the form

$`T_c={\displaystyle \frac{W^2}{6\pi }}{\displaystyle _{\mathrm{}}^\mu }\text{Im}\left[{\displaystyle \frac{g^2}{1\frac{W^2g^2}{6}}}\right]𝑑E`$ (28)
$`={\displaystyle \frac{W\sqrt{2}}{4\pi }}\left[\sqrt{1y^2}{\displaystyle \frac{1}{\sqrt{3}}}\mathrm{tan}^1\sqrt{3(1y^2)}\right],`$ (29)
where $`n=n_0(y)/2`$. If, in addition, $`n1`$, we obtain

$$T_c=\frac{3}{4\sqrt{2}}Wn.53Wx.$$
(30)
Notice, that virtual crystal approximation (for $`J\mathrm{}`$)  gives

$$T_c=\frac{2}{15}Wx.13Wx,$$
(31)
thus being imprecise by a factor of 4.

For arbitrary exchange, the integral in Eq. (24) can be calculated only numerically. The results of such calculation are presented on FIG. 1.

## IV Nonferromagnetism

Note that that our main result (equation for the $`T_c`$) indicates it’s own limits of validity. In part of the $`J/Wn`$ plane, Eq. (24) gives $`T_c<0`$. Negative value of $`T_c`$ means, that at any temperature, including $`T=0`$, the paramagnetic phase is stable with respect to appearance of small spontaneous magnetic moment, and strongly suggests that the ground state in this part of the phase plane is nonferromagnetic . The mechanisms leading to the destruction of ferromagnetism can be easily understood in the extreme cases of weak and strong exchange.

It is well known, that for $`J0`$ the DE Hamiltonian can be reduced to to the RKKY Hamiltonian (see review  and references to the original papers therein)

$`H_{RKKY}={\displaystyle \underset{nn^{}}{}}I_{nn^{}}𝐦_n𝐦_n^{}.`$ (32)
Thus weak exchange between the electrons and the core spins generates effective long-range and sign oscillating effective exchange interaction between the core spins.

In the opposite case $`J\mathrm{}`$, we can reduce the DE Hamiltonian (2) to effective (classic) $`tJ`$ model (see review  and references therein). To do that, first diagonalize the exchange part of the DE Hamiltonian by choosing local spin quantization axis on each site in the direction of $`𝐦`$, the latter being determined by polar angle $`\theta `$ and azimuthal angle $`\varphi `$. In this representation the Hamiltonian is

$`H_{nn^{}}=J\sigma ^z\delta _{nn^{}}+t_{nn^{}}\left(\begin{array}{cc}a_{nn^{}}& b_{nn^{}}\\ b_{n^{}n}^{}& a_{n^{}n}\end{array}\right),`$ (35)
where

$`a_{nn^{}}=\mathrm{cos}{\displaystyle \frac{\theta _n}{2}}\mathrm{cos}{\displaystyle \frac{\theta _n^{}}{2}}+\mathrm{sin}{\displaystyle \frac{\theta _n}{2}}\mathrm{sin}{\displaystyle \frac{\theta _n^{}}{2}}e^{i(\varphi _n^{}\varphi _n)}`$ (36)
$`b_{nn^{}}=\mathrm{sin}{\displaystyle \frac{\theta _n}{2}}\mathrm{cos}{\displaystyle \frac{\theta _n^{}}{2}}e^{i\varphi _n}\mathrm{cos}{\displaystyle \frac{\theta _n}{2}}\mathrm{sin}{\displaystyle \frac{\theta _n^{}}{2}}e^{i\varphi _n^{}}.`$ (37)
The next step is to apply a canonical transformation

$`H\stackrel{~}{H}=e^SHe^S=H+[S,H]+{\displaystyle \frac{1}{2}}[S[S,H]]+\mathrm{}`$ (38)
which excludes all band-to-band transitions. This can be achieved if we chose the operator $`S`$ in the form

$`S_{nn^{}}={\displaystyle \frac{t_{nn^{}}}{2J}}\left(\begin{array}{cc}a_{nn^{}}& b_{nn^{}}\\ b_{n^{}n}^{}& a_{n^{}n}\end{array}\right).`$ (41)
We have

$`[S,H^{ex}]_{nn^{}}=t_{nn^{}}\left(\begin{array}{cc}0& b_{nn^{}}\\ b_{n^{}n}^{}& 0\end{array}\right)`$ (44)
$`[S,[S,H^{ex}]]_{nn}=2{\displaystyle \underset{n^{\prime \prime }}{}}I_{nn^{\prime \prime }}\left(\begin{array}{cc}|b_{nn^{\prime \prime }}|^2& 0\\ 0& |b_{nn^{\prime \prime }}|^2\end{array}\right),`$ (47)
where $`I_{nn^{\prime \prime }}=|t_{nn^{\prime \prime }}|^2/(2J)`$. Keeping terms up to the second order with respect to $`t`$ (and only site-diagonal part of the second order terms) we obtain

$`\stackrel{~}{H}_{nn^{}}=J\sigma ^z\delta _{nn^{}}+t_{nn^{}}\left(\begin{array}{cc}a_{nn^{}}& 0\\ 0& a_{n^{}n}\end{array}\right)`$ (50)
$`+{\displaystyle \underset{n^{\prime \prime }}{}}I_{nn^{\prime \prime }}(𝐦_n𝐦_{n^{\prime \prime }}1)\sigma ^z\delta _{nn^{}},`$ (51)
In the second quantization the Hamiltonian obtained has the form (ignoring the constant term)

$`\stackrel{~}{H}={\displaystyle \underset{nn^{}}{}}t_{nn^{}}[\mathrm{cos}{\displaystyle \frac{\theta _n}{2}}\mathrm{cos}{\displaystyle \frac{\theta _n^{}}{2}}`$ (52)
$`+\mathrm{sin}{\displaystyle \frac{\theta _n}{2}}\mathrm{sin}{\displaystyle \frac{\theta _n^{}}{2}}e^{i(\varphi _n^{}\varphi _n)}]d_n^{}d_n^{}`$ (53)
$`+{\displaystyle \underset{nn^{}}{}}I_{nn^{}}𝐦_n𝐦_n^{}d_n^{}d_n,`$ (54)
where $`d^{}(d)`$ is the operator of creation (annihilation) of a spinless electron. Thus strong exchange between the electrons and the core spins generates effective short-range antiferromagnetic exchange between the core spins.

In conclusion, we formulated the Dynamical Mean Field Approximation equations for the double-exchange model with classical spins for arbitrary electron dispersion law, concentration and relation between the exchange and the electron bandwidth. In the vicinity of the paramagnetic-ferromagnetic transition critical temperature $`T_c`$, these equations were reduced to a MF equation, describing a single spin in an effective field, proportional to the macroscopic magnetization, and we obtained explicit formula for the $`T_c`$. For the semi-circular electron density of states we explicitly calculated the transition temperature as a function of the exchange interaction and electron concentration.

This research was supported by the Israeli Science Foundation administered by the Israel Academy of Sciences and Humanities.


text:
# QCD Effective Actions at High Temperature and the Quasiparticle Picture

## 1 Introduction

The study of the high temperature phase of QCD is necessary to characterize the properties of the quark-gluon plasma (QGP), relevant for the phenomenology of heavy ion collisions. A serious obstacle to make robust predictions of its properties is the non-perturbative character of the long distance physics. The standard Monte Carlo techniques can only be used to get the static properties of the system, either by performing 4d lattice or (dimensional reduced) 3d lattice simulations. Unfortunately, a similar powerful approach to study the dynamics in this phase is still missing. Efforts to develop non-perturbative tools for the QGP are certainly required.

One approach to study dynamical properties of non-Abelian plasmas is based on the fact that the low energy modes behave classically. It is possible to construct effective classical field theories which are amenable to numerical treatment. These theories treat the hard modes as quasiparticles moving in the background of soft classical fields, and they are modeled with a simple transport equation. Previous sessions of SEWM were intensively devoted to the construction of this program, and I refer to their proceedings for the basic references. These effective models allowed one to evaluate transport coefficients to leading order in $`g`$, the gauge coupling constant.

It seems reasonable to ask at which point the predictions of kinetic theory depart from those of quantum field theory. If not, it would be interesting to see whether the analysis explained above could be pushed beyond leading order in $`g`$. With this aim I will review the predictions of kinetic theory for systems close to thermal equilibrium, comparing them with those from quantum field theory. This approach is constructive, as the intention is not to derive transport equations from quantum field theory (while this is possible!), but rather to see if a simple quasiparticle model describes correctly the QCD predictions for the long distance physics. This is based on work done in collaboration with Mikko Laine and Stanisław Mrówczyński. A similar effort in that direction was first carried out in Ref..

There are two different kinetic equations to describe classical colored particles. One of them describes color as a classical degree of freedom. The other one describes color as a quantum degree of freedom, represented by a matrix in a certain representation $`R`$ of the gauge group $`SU(N_c)`$. While the two different kinetic equations predict the same low order operators in these effective theories, they do not for higher order operators.

## 2 Transport Equations for Classical Color

Let us a consider a non-relativistic particle of mass $`m`$ carrying a classical color charge $`Q^a`$, where $`a=1`$ to $`N_c1`$ for a $`SU(N_c)`$ group. The Hamiltonian for this particle is

$$H=\frac{1}{2m}\left(𝐤g𝐀^aQ_a\right)^2+gA_0^aQ_a,$$
(1)
where $`A_\mu ^a`$ is the non-Abelian vector gauge field. The color charge is a dynamical variable, which evolves in time $`t`$, as the canonical momentum $`𝐤`$ and position $`𝐱`$ of the particle. Its phase-space is then enlarged by adding the color charge variables. The Poisson brackets (PB) of this classical system are also modified. Hamiltonian equations of motion can be derived with the help of these Poisson brackets, and they correspond to the well-known Wong equations.

A statistical description of a system of classical colored particles starts with the definition of the one-particle distribution function $`f(𝐱,𝐤,Q;t)`$. In the absence of collisions it evolves according to the Hamiltonian dynamics, $`\frac{df}{dt}=\{H,f\}_{PB}`$. A formulation of the transport equation in terms of the canonical momentum (the phase-space variable) yields a gauge dependent equation222The same problem reappears in quantum field theory derivations of transport equations!. This problem can be overcome by rewriting the equation in terms of the (gauge invariant) kinetic momentum $`𝐩=𝐤g𝐀^aQ_a`$, as this is the variable associated to the particle velocity $`𝐯=𝐩/m`$. Then the equation reads

$$v^\mu \left(D_\mu gQ_aF_{\mu \nu }^a\frac{}{p^\nu }\right)f(x,p,Q)=0,$$
(2)
where $`v^\mu =(1,𝐯)`$, $`D_\mu =_\mu gf^{abc}A_\mu ^bQ^c/Q^a`$, and $`p_0=𝐩^2/2m`$ for a non-relativistic particle.

A Hamiltonian formulation of the dynamics is very useful. In particular, it allows one to recognize constants of motion, and thus getting exact solutions of the collisionless transport equation. If the system is invariant under space translations in the direction $`𝐧`$, then it is easy to prove that $`𝐧𝐤`$ is a constant of motion. While this is not a gauge invariant quantity, it is in the type of gauges that are respectful with the translation invariance of the system, $`n^i_iA_\mu ^a=0`$. Thus, an exact solution to (2) is written as $`f(𝐧(𝐩+g𝐀_aQ^a))`$, for an arbitrary function $`f`$, which should be fixed with a boundary condition.

A generalization of the approach for relativistic systems is straightforward. All what changes in the final equations is the relativistic, rather than non-relativistic, dispersion relation. Using this philosophy exact solutions to the relativistic transport equations in static systems were found in Ref.. The imposed boundary condition was that $`f`$ should reduce to the Fermi-Dirac or Bose-Einstein equilibrium functions when $`g0`$, which then fixes the form of $`f`$. One can also prove that these distribution functions describe equilibrium solutions in the presence of background fields (that is, with no entropy production). The colored current in the plasma is written in terms of the external fields, and one can then compute the associated effective action. When compared to the effective action of a system in the background of static fields as computed in quantum field theory, one notices that it matches for the low dimensional operators (those proportional to the quadratic and cubic Casimirs), but there is a numerical discrepancy in the coefficient of the dimension 4 operator, unless the particles are in a high dimensional color representation. It is possible to identify the reason of this discrepancy in the classical color algebra, and the fact that the $`Q`$ charges are commuting objects. Since quarks and gluons are in low dimensional representations, a quantum treatment of color then seems mandatory for QCD, at least to reach agreement between the two theories for all the static operators.

## 3 Transport Equations for Quantum Color

A (first) quantized treatment of the system described in (1) amounts to replace c-numbers by operators in a Hilbert space. For a particle with color charge in a representation $`R`$ of $`SU(N_c)`$, the color charge $`Q_a`$ is replaced by $`T_a`$, where $`T_a`$ is a matrix, a generator of $`SU(N_c)`$ in the representation $`R`$. Poisson brackets are replaced by commutators. In a Heisenberg representation, one can define the operators of position, momentum and color charge, and derive their equations of motion, according to the Hamiltonian dynamics. These equations reduce to the classical Wong equations in the limit where all operators commute.

The statistical treatment of the system starts with the Wigner operator, defined as a certain Fourier transform of the density matrix. The Wigner operator is also a matrix in color space, whose dimension depends on the representation $`R`$. In a collisionless situation, its dynamical evolution is governed by the Hamiltonian. In the limit where spatial derivatives of higher order can be neglected (that is, performing a gradient expansion), it reduces to

$$[v^\mu D_\mu ,W(p,x)]\frac{g}{2}v^\mu \{F_{\mu \nu }(x),\frac{W(p,x)}{p_\nu }\}=0,$$
(3)
where $`D_\mu =_\mu igA_\mu ^aT^a`$. The equation is the same for non-relativistic and relativistic systems, all what changes is the dispersion relation. These transport equations were first derived from QCD in a gradient expansion in Refs. , while here they have been obtained from a simple quasiparticle model.

In Ref.  exact solutions to (3) were found for systems with a translation invariance in the direction $`n^\mu `$. One could expect that the form of the solutions to be $`W(p,x)=f(n(p+gA^aT^a))`$. However, this is only the case when $`[D(nA),nA]=0`$. If not, gradient (and non-local) terms correct this expression. For static systems, this implies that an equilibrium solution in the presence of background fields is only possible if those are in the direction of the Cartan subalgebra of $`SU(N_c)`$.

One can compare the results arising from a quasiparticle model, made up of quarks, antiquarks and gluons, and those of QCD. One finds a perfect matching in the static limit with constant background fields. A discrepancy is found when compared with dimensional reduced effective theories, if the vector gauge fields do not commute. There is also a subtlety concerning the gluons. The transport equation for the hard gluons needs an infrared cutoff, so as to avoid a double counting of the soft modes (as quasiparticles, as soft classical fields). This cutoff is also needed in order to avoid that the distribution function of the hard gluons becomes negative, loosing its probabilistic interpretation. A matching procedure with the soft classical field modes allows one to eliminate the cutoff dependence in all final physical quantities.

## 4 Conclusions

We have pointed out where the classical and quantum color quasiparticle models agree and disagree. While we have not discussed the structure of the collision terms, a similar discrepancy beyond leading order in $`g`$ might be expected there.

I would like to emphasize that it is worthwhile to explore these quasiparticle models, not only for the reasons mentioned in the introduction, but also because they might teach us a lot about fluctuations around global and local equilibrium in a quark-gluon plasma.


text:
# Instanton model of light hadrons: from low to high energies Talk given at XXXII International Simposium Multiparticle Dynamics, September 7-13, Alushta, Crimea, Ukraine This work is supported by RFBR Grants 01-02-16431 and 02-02-16194 and by Grant INTAS-2000-366.

## Acknowledgments

I am grateful to A.M. Bakulev, W. Broniowski, A. Di Giacomo, A.S. Gorski, N.I. Kochelev, S.V. Mikhailov, M.K. Volkov, L. Tomio and V.L. Yudichev for many useful discussions on topics related to this talk.


text:
# Principal Manifolds and Nonlinear Dimension Reduction via Local Tangent Space Alignment

## 1 Introduction

Many high-dimensional data in real-world applications can be modeled as data points lying close to a low-dimensional nonlinear manifold. Discovering the structure of the manifold from a set of data points sampled from the manifold possibly with noise represents a very challenging unsupervised learning problem . The discovered low-dimensional structures can be further used for classification, clustering, outlier detection and data visualization. Example low-dimensional manifolds embedded in high-dimensional input spaces include image vectors representing the same 3D objects under different camera views and lighting conditions, a set of document vectors in a text corpus dealing with a specific topic, and a set of 0-1 vectors encoding the test results on a set of multiple choice questions for a group of students . The key observation is that the dimensions of the embedding spaces can be very high (e.g., the number of pixels for each images in the image collection, the number of terms (words and/or phrases) in the vocabulary of the text corpus, or the number of multiple choice questions in the test), the intrinsic dimensionality of the data points, however, are rather limited due to factors such as physical constraints and linguistic correlations. Traditional dimension reduction techniques such as principal component analysis and factor analysis usually work well when the data points lie close to a linear (affine) subspace in the input space . They can not, in general, discover nonlinear structures embedded in the set of data points.

Recently, there have been much renewed interests in developing efficient algorithms for constructing nonlinear low-dimensional manifolds from sample data points in high-dimensional spaces, emphasizing simple algorithmic implementation and avoiding optimization problems prone to local minima . Two lines of research of manifold learning and nonlinear dimension reduction have emerged: one is exemplified by  where pairwise geodesic distances of the data points with respect to the underlying manifold are estimated, and the classical multi-dimensional scaling is used to project the data points into a low-dimensional space that best preserves the geodesic distances. Another line of research follows the long tradition starting with self-organizing maps (SOM) , principal curves/surfaces  and topology-preserving networks . The key idea is that the information about the global structure of a nonlinear manifold can be obtained from a careful analysis of the interactions of the overlapping local structures. In particular, the local linear embedding (LLE) method constructs a local geometric structure that is invariant to translations and orthogonal transformations in a neighborhood of each data points and seeks to project the data points into a low-dimensional space that best preserves those local geometries .

Our approach draws inspiration from and improves upon the work in  which opens up new directions in nonlinear manifold learning with many fundamental problems requiring to be further investigated. Our starting point is not to consider nonlinear dimension reduction in isolation as merely constructing a nonlinear projection, but rather to combine it with the process of reconstruction of the nonlinear manifold, and we argue that the two processes interact with each other in a mutually reinforcing way. In this paper, we address two inter-related objectives of nonlinear structure finding: 1) to construct the so-called principal manifold  that goes through “the middle” of the data points; and 2) to find the global coordinate system (the natural parametrization space) that characterizes the set of data points in a low-dimensional space. The basic idea of our approach is to use the tangent space in the neighborhood of a data point to represent the local geometry, and then align those local tangent spaces to construct the global coordinate system for the nonlinear manifold.

The rest of the paper is organized as follows: in section 2, we formulate the problem of manifold learning and dimension reduction in more precise terms, and illustrate the intricacy of the problem using the linear case as an example. In section 3, we discuss the issue of learning local geometry using tangent spaces, and in section 4 we show how to align those local tangent spaces in order to learn the global coordinate system of the underlying manifold. Section 5 discusses how to construct the manifold once the global coordinate system is available. We call the new algorithm local tangent space alignment (LTSA) algorithm. In section 6, we present an error analysis of LTSA, especially illustrating the interactions among curvature information embedded in the Hessian matrices, local sampling density and noise level, and the regularity of the Jacobi matrix. In section 7, we show how the partial eigendecomposition used in global coordinate construction can be efficiently computed. We then present a collection of numerical experiments in section 8. Section 9 concludes the paper and addresses several theoretical and algorithmic issues for further research and improvements.

## 2 Manifold Learning and Dimension Reduction

We assume that a $`d`$-dimensional manifold $``$ embedded in an $`m`$-dimensional space $`(d<m)`$ can be represented by a function

$$f:C^dR^m,$$

where $`C`$ is a compact subset of $`^d`$ with open interior. We are given a set of data points $`x_1,\mathrm{},x_N`$, where $`x_i^m`$ are sampled possibly with noise from the manifold, i.e.,

$$x_i=f(\tau _i)+ϵ_i,i=1,\mathrm{},N,$$

where $`ϵ_i`$ represents noise. By dimension reduction we mean the estimation of the unknown lower dimensional feature vectors $`\tau _i`$’s from the $`x_i`$’s, i.e., the $`x_i`$’s which are data points in $`^m`$ is (nonlinearly) projected to $`\tau _i`$’s which are points in $`^d`$, with $`d<m`$ we realize the objective of dimensionality reduction of the data points. By manifold learning we mean the reconstruction of $`f`$ from the $`x_i`$’s, i.e., for an arbitrary test point $`\tau C^d`$, we can provide an estimate of $`f(\tau )`$. These two problems are inter-related, and the solution of one leads to the solution of the other. In some situations, dimension reduction can be the means to an end by itself, and it is not necessary to learn the manifold. In this paper, however, we promote the notion that both problems are really the two sides of the same coin, and the best approach is not to consider each in isolation. Before we tackle the algorithmic details, we first want to point out that the key difficulty in manifold learning and nonlinear dimension reduction from a sample of data points is that the data points are unorganized, i.e., no adjacency relationship among them are known beforehand. Otherwise, the learning problem becomes the well-researched nonlinear regression problem (for a more detailed discussion, see  where techniques from computational geometry was used to solve error-free manifold learning problems). To ease discussion, in what follows we will call the space where the data points live the input space, and the space into which the data points are projected the feature space.

To illustrate the concepts and problems we have introduced, we consider the example of linear manifold learning and linear dimension reduction. We assume that the set of data points are sampled from a $`d`$-dimensional affine subspace, i.e.,

$$x_i=c+U\tau _i+ϵ_i,i=1,\mathrm{},N,$$

where $`c^m,\tau _i^d`$ and $`ϵ_i^m`$ represents noise. $`U^{m\times d}`$ is a matrix forms an orthonormal basis of the affine subspace. Let

$$X=[x_1,\mathrm{},x_N],T=[\tau _1,\mathrm{},\tau _N],E=[ϵ_1,\mathrm{},ϵ_N].$$

Then in matrix form, the data-generation model can be written as

$$X=ce^T+UT+E,$$

here $`e`$ is an $`N`$-dimensional column vector of all ones. The problem of linear manifold learning amounts to seeking $`c,U`$ and $`T`$ to minimize the reconstruction error $`E`$, i.e.,

$$\mathrm{min}E=\underset{c,U,T}{\mathrm{min}}X(ce^T+UT)_F,$$

where $`_F`$ stands for the Frobenius norm of a matrix. This problem can be readily solved by the singular value decomposition (SVD) based upon the following two observations.

1) The norm of the error matrix $`E`$ can be reduced by removing the mean of the columns of $`E`$ from each column of $`E`$, and hence one can assume that the optimal $`E`$ has zero mean. This requirement can be fulfilled if $`c`$ is chosen as the mean of $`X`$, i.e., $`c=Xe/N\overline{x}`$.

2) The low-rank matrix $`UT`$ is the optimal rank-$`d`$ approximation to the centered data matrix $`X\overline{x}e^T`$. Hence the the optimal solution is given by the SVD of $`X\overline{x}e^T`$,

$$X\overline{x}e^T=Q\mathrm{\Sigma }V^T,P^{m\times m},\mathrm{\Sigma }^{m\times N},V^{N\times N},$$

i.e., $`UT=Q_d\mathrm{\Sigma }_dV_d^T,`$ where $`\mathrm{\Sigma }_d=\mathrm{diag}(\sigma _1,\mathrm{},\sigma _d)`$ with the $`d`$ largest singular values of $`X\overline{x}e^T`$, $`Q_d`$ and $`V_d`$ are the matrices of the corresponding left and right singular vectors, respectively. The optimal $`U^{}`$ is then given by $`Q_d`$ and the learned linear manifold is represented by the linear function

$$f(\tau )=\overline{x}+U^{}\tau .$$

In this model, the coordinate matrix $`T`$ corresponding to the data matrix $`X`$ is given by

$$T=(U^{})^T(X\overline{x}e^T)=\mathrm{diag}(\sigma _1,\mathrm{},\sigma _d)V_d^T.$$

Ideally, the dimension $`d`$ of the learned linear manifold should be chosen such that $`\sigma _{d+1}\sigma _d`$.

The function $`f`$ is not unique in the sense that it can be reparametrized, i.e., the coordinate can be replaced by $`\stackrel{~}{\tau }`$ with a global affine transformation $`\tau =P\stackrel{~}{\tau }`$, if we change the basis matrix $`U^{}`$ to $`U^{}P`$. What we are interested in with respect to dimension reduction is the low-dimensional representation of the linear manifold in the feature space. Therefore, without loss of generality, we can assume that the feature vectors are uniformly distributed. For a given data set, this amounts to assuming that the coordinate matrix $`T`$ is orthonormal in row, i.e., $`TT^T=I_d`$. Hence we we can take $`T=V_d^T`$ and the linear function is now the following

$$f(\tau )=\overline{x}+U^{}\mathrm{diag}(\sigma _1,\mathrm{},\sigma _d)\tau .$$

For the linear case we just discussed, the problem of dimension reduction is solved by computing the right singular vectors $`V_d`$, and this can be done without the help of the linear function $`f`$. Similarly, the construction of the linear function $`f`$ is done by computing $`U^{}`$ which are just the $`d`$ largest left singular vectors of $`X\overline{x}e^T`$.

The case for nonlinear manifolds is more complicated. In general, the global nonlinear structure will have to come from local linear analysis and alignment . In , local linear structure of the data set are extracted by representing each point $`x_i`$ as a weighted linear combination of its neighbors, and the local weight vectors are preserved in the feature space in order to obtain a global coordinate system. In , a linear alignment strategy was proposed for aligning a general set of local linear structures. The type of local geometric information we use is the tangent space at a given point which is constructed from a neighborhood of the given point. The local tangent space provides a low-dimensional linear approximation of the local geometric structure of the nonlinear manifold. What we want to preserve are the local coordinates of the data points in the neighborhood with respect to the tangent space. Those local tangent coordinates will be aligned in the low dimensional space by different local affine transformations to obtain a global coordinate system. Our alignment method is similar in spirit to that proposed in . In the next section we will discuss the local tangent space and global alignment which will then be applied to data points sampled with noise in Section 4.

## 3 Local Tangent Space and Its Global Alignment

We assume that $``$ is a $`d`$-dimensional manifold in a $`m`$-dimensional space with unknown generating function $`f(\tau ),\tau R^d`$, and we are given a data set consists of $`N`$ $`m`$-dimensional vectors $`X=[x_1,\mathrm{},x_N],x_i^m`$ generated from the following noise-free model,

$$x_i=f(\tau _i),i=1,\mathrm{},N,$$

where $`\tau _i^d`$ with $`d<m`$. The objective as we mentioned before for nonlinear dimension reduction is to reconstruct $`\tau _i`$’s from the corresponding function values $`f(\tau _i)`$’s without explicitly constructing $`f`$. Assume that the function $`f`$ is smooth enough, using first-order Taylor expansion at a fixed $`\tau `$, we have

$`f(\overline{\tau })=f(\tau )+J_f(\tau )(\overline{\tau }\tau )+O(\overline{\tau }\tau ^2),`$ (3.1)
where $`J_f(\tau )^{m\times d}`$ is the Jacobi matrix of $`f`$ at $`\tau `$. If we write the $`m`$ components of $`f(\tau )`$ as

$$f(\tau )=\left[\begin{array}{c}f_1(\tau )\\ \mathrm{}\\ f_m(\tau )\end{array}\right],\text{then }J_f(\tau )=\left[\begin{array}{ccc}f_1/\tau _1& \mathrm{}& f_1/\tau _d\\ \mathrm{}& \mathrm{}& \mathrm{}\\ f_m/\tau _1& \mathrm{}& f_m/\tau _d\end{array}\right].$$

The tangent space $`𝒯_\tau `$ of $`f`$ at $`\tau `$ is spanned by the $`d`$ column vectors of $`J_f(\tau )`$ and is therefore of dimension at most $`d`$, i.e., $`𝒯_\tau =\mathrm{span}(J_f(\tau ))`$. The vector $`\tau \overline{\tau }`$ gives the coordinate of $`f(\tau )`$ in the affine subspace $`f(\tau )+𝒯_\tau `$. Without knowing the function $`f`$, we can not explicitly compute the Jacobi matrix $`J_f(\tau )`$. However, if we know $`𝒯_\tau `$ in terms of $`Q_\tau `$, a matrix forming an orthonormal basis of $`𝒯_\tau `$, we can write

$$J_f(\tau )(\overline{\tau }\tau )=Q_\tau \theta _\tau ^{},$$

Furthermore,

$$\theta _\tau ^{}=Q_\tau ^TJ_f(\tau )(\tau \overline{\tau })P_\tau (\overline{\tau }\tau ).$$

The mapping from $`\tau `$ to $`\theta _\tau ^{}`$ represents a local affine transformation. This affine transformation is unknown because we do not know the function $`f`$. The vector $`\theta _\tau ^{}`$, however, has an approximate $`\theta _\tau `$ that orthogonally projects $`f(\overline{\tau })f(\tau )`$ onto $`𝒯_\tau `$,

$$\theta _\tau Q_\tau ^T(f(\overline{\tau })f(\tau ))=\theta _\tau ^{}+O(\overline{\tau }\tau ^2),$$
(3.2)
provided $`Q_\tau `$ is known at each $`\tau `$. Ignoring the second-order term, the global coordinate $`\tau `$ satisfies

$$𝑑\tau _{\mathrm{\Omega }(\tau )}P_\tau (\overline{\tau }\tau )\theta _\tau 𝑑\overline{\tau }0.$$

Here $`\mathrm{\Omega }(\tau )`$ defines the neighborhood of $`\tau `$. Therefore, a natural way to approximate the global coordinate is to find a global coordinate $`\tau `$ and a local affine transformation $`P_\tau `$ that minimize the error function

$`{\displaystyle 𝑑\tau _{\mathrm{\Omega }(\tau )}P_\tau (\overline{\tau }\tau )\theta _\tau 𝑑\overline{\tau }}.`$ (3.3)
This represents a nonlinear alignment approach for the dimension reduction problem (this idea will be picked up at the end of section 4).

On the other hand, a linear alignment approach can be devised as follows. If $`J_f(\tau )`$ is of full column rank, the matrix $`P_\tau `$ should be non-singular and

$$\overline{\tau }\tau P_\tau ^1\theta _\tau L_\tau \theta _\tau .$$

The above equation shows that the affine transformation $`L_\tau `$ should align this local coordinate with the global coordinate $`\tau \overline{\tau }`$ for $`f(\tau )`$. Naturally we should seek to find a global coordinate $`\tau `$ and a local affine transformation $`L_\tau `$ to minimize

$`{\displaystyle 𝑑\tau _{\mathrm{\Omega }(\tau )}\overline{\tau }\tau L_\tau \theta _\tau 𝑑\overline{\tau }}.`$ (3.4)
The above amounts to matching the local geometry in the feature space. Notice that $`\theta _\tau `$ is defined by the “known” function value and the “unknown” orthogonal basis matrix $`Q_\tau `$ of the tangent space. It turns out, however, $`Q_\tau `$ can be approximately determined by certain function values. We will discuss this approach in the next section. Clearly, this linear approach is more readily applicable than (3.3). Obviously, If the manifold $``$ is not regular, i.e., the Jacobi matrix $`J_f`$ is not of full column rank at some points $`\tau C`$, then the two minimization problems (3.4) and (3.3) may lead to quite different solutions.

As is discussed in the linear case, the low-dimensional feature vector $`\tau `$ is not uniquely determined by the manifold $``$. We can reparametrize $``$ using $`f(g(\tau ))`$ where $`g()`$ is a smooth $`1`$-to-$`1`$ onto mapping of $`C`$ to itself. The parameterization of $``$ can be fixed by requiring that $`\tau `$ has a uniform distribution over $`C`$. This will come up as a normalization issue in the next section.

## 4 Feature Extraction through Alignment

Now we consider how to construct the global coordinates and local affine transformation when we are given a data set $`X=[x_1,\mathrm{},x_N]`$ sampled with noise from an underlying nonlinear manifold,

$$x_i=f(\tau _i)+ϵ_i,i=1,\mathrm{},N,$$

where $`\tau _i^d`$, $`x_i^m`$ with $`d<m`$. For each $`x_i`$, let $`X_i=[x_{i_1},\mathrm{},x_{i_k}]`$ be a matrix consisting of its $`k`$-nearest neighbors including $`x_i`$, say in terms of the Euclidean distance. Consider computing the best $`d`$-dimensional affine subspace approximation for the data points in $`X_i`$,

$$\underset{x,\mathrm{\Theta },Q}{\mathrm{min}}\underset{j=1}{\overset{k}{}}x_{i_j}(x+Q\theta _j)_2^2=\underset{x,\mathrm{\Theta },Q}{\mathrm{min}}X_i(xe^T+Q\mathrm{\Theta })_2^2,$$

where $`Q`$ is of $`d`$ columns and is orthonormal, and $`\mathrm{\Theta }=[\theta _1,\mathrm{},\theta _k]`$. As is discussed in section 2, the optimal $`x`$ is given by $`\overline{x}_i`$, the mean of all the $`x_{i_j}`$’s and the optimal $`Q`$ is given by $`Q_i`$, the $`d`$ left singular vectors of $`X_i(Iee^T/k)`$ corresponding to its $`d`$ largest singular values, and $`\mathrm{\Theta }`$ is given by $`\mathrm{\Theta }_i`$ defined as

$$\mathrm{\Theta }_i=Q_i^TX_i(I\frac{1}{k}ee^T)=[\theta _1^{(i)},\mathrm{},\theta _k^{(i)}],\theta _j^{(i)}=Q_i^T(x_{i_j}\overline{x}_i).$$
(4.5)
Therefore we have

$$x_{i_j}=\overline{x}_i+Q_i\theta _j^{(i)}+\xi _j^{(i)},$$
(4.6)
where $`\xi _j^{(i)}=(IQ_iQ_i^T)(x_{i_j}\overline{x}_i)`$ denotes the reconstruction error.

We now consider constructing the global coordinates $`\tau _i,i=1,\mathrm{},N`$, in the low-dimensional feature space based on the local coordinates $`\theta _j^{(i)}`$ which represents the local geometry. Specifically, we want $`\tau _{i_j}`$ to satisfy the following set of equations, i.e., the global coordinates should respect the local geometry determined by the $`\theta _j^{(i)}`$,

$`\tau _{i_j}=\overline{\tau }_i+L_i\theta _j^{(i)}+ϵ_j^{(i)},j=1,\mathrm{},k,i=1,\mathrm{},N,`$ (4.7)
where $`\overline{\tau }_i`$ is the mean of $`\tau _{i_j},j=1,\mathrm{},k`$. In matrix form,

$$T_i=\frac{1}{k}T_iee^T+L_i\mathrm{\Theta }_i+E_i,$$

where $`T_i=[\tau _{i_1},\mathrm{},\tau _{i_k}]`$ and $`E_i=[ϵ_1^{(i)},\mathrm{},ϵ_k^{(i)}]`$ is the local reconstruction error matrix, and we write

$`E_i=T_i(I{\displaystyle \frac{1}{k}}ee^T)L_i\mathrm{\Theta }_i.`$ (4.8)
To preserve as much of the local geometry in the low-dimensional feature space, we seek to find $`\tau _i`$ and the local affine transformations $`L_i`$ to minimize the reconstruction errors $`ϵ_j^{(i)}`$, i.e.,

$$\underset{i}{}E_i^2\underset{i}{}T_i(I\frac{1}{k}ee^T)L_i\mathrm{\Theta }_i^2=\mathrm{min}.$$
(4.9)
Obviously, the optimal alignment matrix $`L_i`$ that minimizes the local reconstruction error $`E_i_F`$ for a fixed $`T_i`$, is given by

$$L_i=T_i(I\frac{1}{k}ee^T)\mathrm{\Theta }_i^+=T_i\mathrm{\Theta }_i^+,\text{and therefore}E_i=T_i(I\frac{1}{k}ee^T)(I\mathrm{\Theta }_i^+\mathrm{\Theta }_i),$$

where $`\mathrm{\Theta }_i^+`$ is the Moor-Penrose generalized inverse of $`\mathrm{\Theta }_i`$. Let $`T=[\tau _1,\mathrm{},\tau _N]`$ and $`S_i`$ be the $`0`$-$`1`$ selection matrix such that $`TS_i=T_i,`$. We then need to find $`T`$ to minimize the overall reconstruction error

$$\underset{i}{}E_i_F^2=TSW_F^2,$$

where $`S=[S_1,\mathrm{},S_N]`$, and $`W=\mathrm{diag}(W_1,\mathrm{},W_N)`$ with

$$W_i=(I\frac{1}{k}ee^T)(I\mathrm{\Theta }_i^+\mathrm{\Theta }_i).$$
(4.10)
To uniquely determine $`T`$, we will impose the constraints $`TT^T=I_d`$, it turns out that the vector $`e`$ of all ones is an eigenvector of

$`BSWW^TS^T`$ (4.11)
corresponding to a zero eigenvalue, therefore, the optimal $`T`$ is given by the $`d`$ eigenvectors of the matrix $`B`$, corresponding to the $`2`$nd to $`d`$+1st smallest eigenvalues of $`B`$.

Remark. We now briefly discuss the nonlinear alignment idea mentioned in (3.3). In particular, in a neighborhood of a data point $`x_i`$ consisting of data points $`X_i=[x_{i_1},\mathrm{},x_{i_k}]`$, by first order Taylor expansion, we have

$$X_i(Iee^T/k)J_f^{(i)}T_i(Iee^T/k).$$

Let $`S_i`$ be the neighborhood selection matrix as defined before, we seek to find $`J_f^{(i)}^{m\times d}`$ and $`T`$ to minimize

$$E(J,T)\underset{i=1}{\overset{N}{}}(XJ_f^{(i)}T)S_i(Iee^T/k)_F^2,$$

where $`J=[J_f^{(1)},\mathrm{},J_f^{(N)}]`$. The LTSA algorithm can be considered as an approach to find an approximate solution to the above minimization problem. We can, however, seek to find the optimal solution of $`E(J,T)`$ using an alternating least squares approach: fix $`J`$ minimize $`E(J,T)`$ with respect to $`T`$, and fix $`T`$ minimize $`E(J,T)`$ with respect to $`J`$, and so on. As an initial value to start the alternating least squares, we can use the $`T`$ obtained from the LTSA algorithm. The details of the algorithm will be presented in a separate paper.

Remark. The minimization problem (4.9) needs certain constraints (i.e., normalization conditions) to be well-posed, otherwise, one can just choose both $`T_i`$ and $`L_i`$ to be zero. However, there are more than one way to impose the normalization conditions. The one we have selected, i.e., $`TT^T=I_d`$, is just one of the possibilities. To illustrate the issue we look at the following minimization problem,

$$\underset{X,Y}{\mathrm{min}}XYA_F$$

The approach we have taken amounts to substituting $`Y=XA^+`$, and minimize $`X(IA^+A)_F`$ with the normalization condition $`XX^T=I`$. However,

$$XYA_F=[X,Y]\left[\begin{array}{c}I\\ A\end{array}\right]_F,$$

and we can minimize the above by imposing the normalization condition

$$[X,Y][X,Y]^T=I.$$

This nonuniqueness issue is closely related to to nonuniqueness of the parametrization of the nonlinear manifold $`f(\tau )`$, which can reparametrized as $`f(\tau (\eta ))`$with a 1-to-1 mapping $`\tau (\eta )`$.

## 5 Constructing Principal Manifolds

Once the global coordinates $`\tau _i`$ are computed for each of the data points $`x_i`$, we can apply some non-parametric regression methods such as local polynomial regression to $`\{(\tau _i,x_i)\}_{i=1}^N`$ to construct the principal manifold underlying the set of points $`x_i`$. Here each of the component functions $`f_j(\tau )`$ can be constructed separately, for example, we have used the simple loess function  in some of our experiments for generating the principal manifolds.

In general, when the low-dimensional coordinates $`\tau _i`$ are available, we can construct an mapping from the $`\tau `$-space (feature space) to the $`x`$-space (input space) as follows.

1. For each fixed $`\tau `$, let $`\tau _i`$ be the nearest neighbor (i.e., $`\tau \tau _i\tau \tau _j,`$ for $`ji`$). Define

$$\theta =L_i^1(\tau \overline{\tau }_i),$$

where $`\overline{\tau }_i`$ be the mean of the feature vectors in a neighbor to which $`\tau _i`$ belong.

2. Back in the input space, we define

$$x=\overline{x}_i+Q_i\theta .$$

Let us define by $`g:\tau x`$ the resulted mapping,

$`g(\tau )=\overline{x}_i+Q_iL_i^1(\tau \overline{\tau }_i).`$ (5.12)
To distinguish the computed coordinates $`\tau _i`$ from the exact ones, in the rest of this paper, we denote by $`\tau _i^{}`$ the exact coordinate, i.e.,

$`x_i=f(\tau _i^{})+ϵ_i^{}.`$ (5.13)
Obviously, the errors of the reconstructed manifold represented by $`g`$ depend on the sample errors $`ϵ_i^{}`$, the local tangent subspace reconstruction errors $`\xi _j^{(i)}`$, and the alignment errors $`ϵ_j^{(i)}`$. The following result show that this dependence is linear.

###### Theorem 5.1

Let $`ϵ_i^{}=x_if(\tau _i^{})`$, $`\xi _j^{(i)}=(IQ_iQ_i^T)(x_i\overline{x}_i)`$, and $`ϵ_i=\tau _i\overline{\tau }_iL_iQ_i^T(x_i\overline{x}_i)`$. Then

$$g(\tau _i)f(\tau _i^{})_2ϵ_i^{}_2+\xi _i_2+L_i^1ϵ_i_2.$$

###### Proof 5.2.

Substituting $`L_i^1(\tau _i\overline{\tau }_i)=L_i^1ϵ_i+Q_i^T(x_i\overline{x}_i)`$ into (5.12) gives

$`g(\tau _i)`$ $`=`$ $`\overline{x}_i+Q_iL_i^1(\tau _i\overline{\tau }_i)`$
$`=`$ $`\overline{x}_i+Q_iQ_i^T(x_i\overline{x}_i)+Q_iL_i^1ϵ_1^{(i)}.`$
Because $`Q_iQ_i^T(x_i\overline{x}_i)=x_i\overline{x}_i\xi _j^{(i)}`$, we obtain that

$`g(\tau _i)`$ $`=`$ $`x_i\xi _j^{(i)}+Q_iL_i^1ϵ_1^{(i)}`$
$`=`$ $`f(\tau _i^{})+ϵ_i^{}\xi _j^{(i)}+Q_iL_i^1ϵ_1^{(i)}.`$
Therefore we have

$$g(\tau _i)f(\tau _i^{})_2ϵ_i^{}_2+\xi _i_2+L_i^1ϵ_i_2,$$

completing the proof.

In the next section, we will give a detailed error analysis to estimate the errors of alignment and tangent space approximation in terms of the noise, the geometric properties of the generating function $`f`$ and the density of the generating coordinates $`\tau _i^{}`$.

## 6 Error Analysis

As is mentioned in the previous section, we assume that that the data points are generated by

$$x_i=f(\tau _i^{})+ϵ_i^{},i=1,\mathrm{},N.$$

For each $`x_i`$, let $`X_i=[x_{i_1},\mathrm{},x_{i_k}]`$ be a matrix consisting of its $`k`$-nearest neighbors including $`x_i`$ in terms of the Euclidean distance. Similar to $`E_i`$ defined in (4.8), we denote by $`E_i^{}`$ the corresponding local noise matrix, $`E_i^{}=[ϵ_{i_1}^{},\mathrm{},ϵ_{i_k}^{}]`$. The low-dimensional embedding coordinate matrix computed by the LTSA algorithm is denoted by $`T=[\tau _1,\mathrm{},\tau _N]`$. We first present a result that bounds $`E_i`$ in terms of $`E_i^{}.`$

###### Theorem 6.1.

Assume $`T^{}=[\tau _1^{},\mathrm{},\tau _N^{}]`$ satisfies $`(T^{})^TT^{}=U_d`$. Let $`\overline{\tau }_i`$ be the mean of $`\tau _{i_1},\mathrm{},\tau _{i_k}`$, Denote $`P_i=Q_i^TJ_f(\overline{\tau }_i^{})`$ and $`H_f_{\mathrm{}}(\overline{\tau }_i^{})`$ the Hessian matrix of the $`\mathrm{}`$-th component function of $`f`$. If the $`P_i`$’s are nonsingular, then

$$E_i_FP_i^1_F(\delta _i+E_i^{}_F),$$

where $`\delta _i`$ is defined by

$$\delta _i^2=\underset{\mathrm{}=1}{\overset{m}{}}\underset{j=1}{\overset{k}{}}H_f_{\mathrm{}}(\overline{\tau }_i^{})_2^2\tau _{i_j}^{}\overline{\tau }_i^{}_2^4$$

Furthermore, if each neighborhood is of size $`O(\eta )`$, then $`EP_i^1_FE^{}+O(\eta ^2)`$.

###### Proof 6.2.

First by (4.8), we have

$`E_i=T_i(I{\displaystyle \frac{1}{k}}ee^T)L_i\mathrm{\Theta }_i=(T_iL_iQ_i^TX_i)(I{\displaystyle \frac{1}{k}}ee^T).`$ (6.14)
To represent $`X_i`$ in terms of the Jacobi matrix of $`f`$, we assume that $`f`$ is smooth enough and use Taylor expansion at $`\overline{\tau }_i^{}`$, the mean of the $`k`$ neighbors of $`\tau _i^{}`$, we have

$$x_{i_j}=f(\overline{\tau }_i^{})+J_i(\tau _{i_j}^{}\overline{\tau }_i^{})+\delta _j^{(i)}+ϵ_{i_j},$$

where $`J_i=J_f(\overline{\tau }_i^{})`$ and $`\delta _j^{(i)}`$ represents the remainder term beyond the first order expansion, in particular, its $`\mathrm{}`$-th components can be approximately written as (using second order approximation),

$$\delta _{\mathrm{},j}^{(i)}\frac{1}{2}(\tau _{i_j}^{}\overline{\tau }_i^{})^TH_f_{\mathrm{}}(\overline{\tau }_i^{})(\tau _{i_j}^{}\overline{\tau }_i^{})$$

with the Hessian matrix $`H_f_{\mathrm{}}(\overline{\tau }_i^{})`$ of the $`\mathrm{}`$-th component function $`f_{\mathrm{}}`$ of $`f`$ at $`\overline{\tau }_i^{}`$. We have in matrix form,

$$X_i=f(\overline{\tau }_i^{})e^T+J_iT_i^{}(I\frac{1}{k}ee^T)+\mathrm{\Delta }_i+E_i^{}$$

with $`\mathrm{\Delta }_i=[\delta _1^{(i)},\mathrm{},\delta _k^{(i)}]`$. Multiplying by the centering matrix $`I\frac{1}{k}ee^T`$ gives

$`X_i(I{\displaystyle \frac{1}{k}}ee^T)=(J_iT_i^{}+\mathrm{\Delta }_i+E_i^{})(I{\displaystyle \frac{1}{k}}ee^T).`$ (6.15)
Substituting (6.15) into (6.14) and denoting $`P_i=Q_i^TJ_i`$, we obtain that

$`E_i=(T_iL_iP_iT_i^{}L_iQ_i^T(\mathrm{\Delta }_i+E_i^{}))(I{\displaystyle \frac{1}{k}}ee^T).`$ (6.16)
For any $`\stackrel{~}{T}`$ satisfying the orthogonal condition $`\stackrel{~}{T}\stackrel{~}{T}^T=I_d`$ and any $`\stackrel{~}{L}_i`$, we also have the similar expression of (6.16) for $`\stackrel{~}{T}_i`$ and $`\stackrel{~}{L}_i`$. Note that $`T`$ and $`L_i`$, $`i=1,\mathrm{},N`$, minimize the overall reconstruction error, $`E_F\stackrel{~}{E}_F`$. Setting $`\stackrel{~}{T}=T^{}`$ and $`\stackrel{~}{L}_i=P_i^1`$, we obtain the upper bound

$$E_i_FP_i^1_2(\mathrm{\Delta }_i_F+E_i^{}_F).$$

We estimate the norm $`\mathrm{\Delta }_i_F`$ by ignoring the higher order terms, and obtain that

$$\mathrm{\Delta }_i_F^2\underset{\mathrm{}=1}{\overset{m}{}}\underset{j=1}{\overset{k}{}}H_f_{\mathrm{}}(\overline{\tau }_i^{})_2^2\tau _{i_j}^{}\overline{\tau }_i^{}_2^4=\delta ^2,$$

completing the proof.

The non-singularity of the matrix $`P_i`$ requires that the Jacobi matrix $`J_i`$ be of full column rank and the two subspaces $`\mathrm{span}(J_i)`$ and the $`d`$ largest left singular vector space $`\mathrm{span}(Q_i)`$ are not orthogonal to each other. We now give a quantitative measurement of the non-singularity of $`P_i`$.

###### Theorem 6.3.

Let $`\sigma _d(\stackrel{~}{J}_i)`$ be the $`d`$-th singular value of $`\stackrel{~}{J}_iJ_iT_i^{}(I\frac{1}{k}ee^T)`$, and denote $`\alpha _i=4(E_i^{}_F+\delta _i)/\sigma _d(\stackrel{~}{J}_i)`$ with $`\delta _i`$ defined in Theorem 6.1. Then

$$P_i^1_F(1+\alpha _i^2)^{1/2}J_i^+_F.$$

###### Proof 6.4.

The proof is simple. Let $`\stackrel{~}{J}_i=U_J\mathrm{\Sigma }_JV_J^T`$ be the SVD of the matrix $`\stackrel{~}{J}_i`$. By (6.15) and perturbation bounds for singular subspaces \[5, Theorem 8.6.5\], the singular vector matrix $`Q_i`$ can be expressed as

$`Q_i=(U_J+U_J^{}H)(I+H_i^TH_i)^{1/2}`$ (6.17)
with

$$H_i_F\frac{4}{\sigma _d(\stackrel{~}{J}_i)}(E^{}_F+\mathrm{\Delta }_i_F))\alpha _i,$$

where $`\sigma _d(\stackrel{~}{J}_i)`$ is the $`d`$-largest singular value of $`\stackrel{~}{J}_i`$. On the other hand, from the SVD of $`\stackrel{~}{J}_i`$, we have $`J_iT_i^{}V_J=U_J\mathrm{\Sigma }_J`$, which gives

$$J_i=U_J\mathrm{\Sigma }_J\left(T_i^{}V_J\right)^1.$$

It follows that

$$P_i=Q_i^TJ_i=(I+H_i^TH_i)^{1/2}\mathrm{\Sigma }_J\left(T_i^{}V_J\right)^1=(I+H_i^TH_i)^{1/2}U_J^TJ_i.$$

Therefore we have

$$P_i^1_F(1+H_i_F^2)^{1/2}J_i^+_F,$$

completing the proof.

The degree of non-singularity of $`J_i`$ is determined by the curvature of the manifold and the rotation of the singular subspace is mainly affected by the sample noises $`ϵ_j`$’s and the neighborhood structure of $`x_i`$’s. The above error bounds clearly show that reconstruction accuracy will suffer if the manifold underlying the data set has singular or near-singular points. This phenomenon will be illustrated in the numerical examples in section 8. Finally, we give an error upper bound for the tangent subspace approximation.

###### Theorem 6.5.

Let $`\mathrm{cond}(\stackrel{~}{J}_i)=\sigma _1(\stackrel{~}{J}_i)/\sigma _d(\stackrel{~}{J}_i)`$ be the spectrum condition number of the $`d`$-column matrix $`\stackrel{~}{J}_i`$. Then

$$X_i(\overline{x}_ie^T+Q_i\mathrm{\Theta }_i)_F\left(1+4(1+\alpha _i^2)\mathrm{cond}(\stackrel{~}{J}_i)\right)(E_i^{}_F+\delta _i).$$

###### Proof 6.6.

By (6.15), we write

$$(IQ_iQ_i^T)X(I\frac{1}{k}ee^T)=(IQ_iQ_i^T)\stackrel{~}{J}_i+\stackrel{~}{\mathrm{\Delta }}_i,$$

with $`\stackrel{~}{\mathrm{\Delta }}_i_FE_i^{}_F+\delta _i`$. To estimate $`(IQ_iQ_i^T)\stackrel{~}{J}_i_F`$, we use the expression (6.17) to obtain

$`(IQ_iQ_i^T)\stackrel{~}{J}_i`$ $`=`$ $`U\left(\left(\begin{array}{cc}I& \\ O& \end{array}\right)\left(\begin{array}{cc}I& \\ H_i& \end{array}\right)(I+H_i^TH_i)^1\right)\mathrm{\Sigma }_JV_J^T`$
$`=`$ $`U\left(\begin{array}{cc}H_i^T& \\ I& \end{array}\right)H_i(I+H_i^TH_i)^1\mathrm{\Sigma }_JV_J^T.`$
Taking norms gives that

$$(IQ_iQ_i^T)\stackrel{~}{J}_i_F(1+H_i_2^2)H_i_F\stackrel{~}{J}_i_24(1+\alpha _i^2)(E_i^{}_F+\delta _i)\mathrm{cond}(\stackrel{~}{J}_i).$$

The result required follows.

The above results show that the accurate determination of the local tangent space is dependent on several factors: curvature information embedded in the Hessian matrices, local sampling density and noise level, and the regularity of the Jacobi matrix.

## 7 Numerical Computation Issues

One major computational cost of LTSA involves the computation of the smallest eigenvectors of the symmetric positive semi-defined matrix $`B`$ defined in (4.11). $`B`$ in general will be quite sparse because of the local nature of the construction of the neighborhoods. Algorithms for computing a subset of the eigenvectors for large and/or sparse matrices are based on computing projections of $`B`$ onto a sequence of Krylov subspaces of the form

$$K_p(B,v_0)=\mathrm{span}\{v_0,Bv_0,B^2v_0,\mathrm{},B^{p1}v_0\},$$

for some initial vectors $`v_0`$ . Hence the computation of matrix-vector multiplications $`Bx`$ needs to be done efficiently. Because of the special nature of $`B`$, $`Bx`$ can be computed neighborhood by neighborhood without explicitly forming $`B`$,

$$Bx=S_1W_1W_1^TS_1^Tx+\mathrm{}+S_NW_NW_N^TS_N^Tx,$$

where as defined in (4.10),

$$W_i=(I\frac{1}{k})(I\mathrm{\Theta }_i^+\mathrm{\Theta }_i).$$

Each term in the above summation only involves the $`x_i`$’s in one neighborhood.

The matrix $`\mathrm{\Theta }_i^+\mathrm{\Theta }_i`$ in the right factor of $`W_i`$ is the orthogonal projector onto the subspace spanned by the rows of $`(\mathrm{\Theta }_i)`$. If the SVD of $`X_i\overline{x}_ie^T=Q^{(i)}\mathrm{\Sigma }^{(i)}(H^{(i)})^T`$ is available, the orthogonal projector is given by $`\mathrm{\Theta }_i^+\mathrm{\Theta }_i=H_iH_i^T`$, where $`H_i`$ is the submatrix of first $`d`$ columns of $`H^{(i)}`$. Otherwise, one can compute the QR decomposition $`\mathrm{\Theta }_i^T=H_iR_i`$ of $`\mathrm{\Theta }_i^T`$ and obtain $`\mathrm{\Theta }_i^+\mathrm{\Theta }_i=H_iH_i^T`$ . Clearly, we have $`H^Te=0`$ because $`\mathrm{\Theta }_ie=0`$. Then we can rewrite $`W_i`$ as

$$W_i=I\frac{1}{k}ee^TH_iH_i^T=I[e/\sqrt{k},H_i][e/\sqrt{k},H_i]^TIG_iG_i^T.$$

It is a orthogonal projector onto the null space spanned by the rows of $`\mathrm{\Theta }_i`$ and $`e^T/\sqrt{k}`$. Therefore the matrix-vector product $`y=S_iW_iW_i^TS_i^Tx`$ can be easily computed as follows: denote by $`I_i=\{i_1,\mathrm{},i_k\}`$ the set of indices for the $`k`$ nearest-neighbors of $`x_i`$, then the $`j`$-th element of $`y`$ is zero for $`jI_i`$ and

$$y(I_i)=x(I_i)G_i(G_i^Tx(I_i)).$$

Here $`y(I_i)=[y_{i_1},\mathrm{},y_{i_k}]^T`$ denotes the section of $`y`$ determined by the neighborhood set $`I_i`$.

If one needs to compute the $`d`$ smallest eigenvectors that are orthogonal to $`e`$ by applying some eigen-solver with an explicitly formed $`B`$. The matrix $`B`$ can be computed by carrying out a partial local summation as follows

$`B(I_i,I_i)B(I_i,I_i)+IG_iG_i^T,i=1,\mathrm{},N`$ (7.20)
with initial $`B=0`$.

Now we are ready to present our Local Tangent Space Alignment (LTSA) algorithm.

Algorithm LTSA (Local Tangent Space Alignment). Given $`N`$ $`m`$-dimensional points sampled possibly with noise from an underlying $`d`$-dimensional manifold, this algorithm produces $`N`$ $`d`$-dimensional coordinates $`T^{d\times N}`$ for the manifold constructed from $`k`$ local nearest neighbors.    \[Extracting local information.\] For each $`i=1,\mathrm{},N`$,    Determine $`k`$ nearest neighbors $`x_{i_j}`$ of $`x_i,j=1,\mathrm{},k`$.    Compute the $`d`$ largest eigenvectors $`g_1,\mathrm{},g_d`$ of the correlation matrix $`(X_i\overline{x}_ie^T)^T(X_i\overline{x}_ie^T)`$, and set
$$G_i=[e/\sqrt{k},g_1,\mathrm{},g_d].$$
\[Constructing alingment matrix.\] Form the matrix $`B`$ by locally summing (7.20) if a direct eigen-solver will be used. Otherwise implement a routine that computes matrix-vector multiplication $`Bu`$ for an arbitrary vector $`u`$.    \[Aligning global cordinates.\] Compute the $`d+1`$ smallest eigenvectors of $`B`$ and pick up the eigenvector matrix $`[u_2,\mathrm{},u_{d+1}]`$ corresponding to the 2nd to $`d+1`$st smallest eigenvalues, and set $`T=[u_2,\mathrm{},u_{d+1}]^T`$.

## 8 Experimental Results

In this section, we present several numerical examples to illustrate the performance of our LTSA algorithm. The test data sets include curves in 2D/3D Euclidean spaces, and surfaces in 3D Euclidean spaces. Especially, we take a closer look at the effects of singular points of a manifold and the interaction of noise levels and sample density. To show that our algorithm can also handle data points in high-dimensional spaces, we also consider curves and surfaces in Euclidean spaces with dimension equal to $`100`$ and an image data set with dimension $`4096`$.

First we test our LTSA method for 1D manifolds (curves) in both 2D and 3D. For a given 1D manifold $`f(\tau )`$ with uniformly sampled coordinates $`\tau _1^{},\mathrm{},\tau _N^{}`$ in a fixed interval, we add Gaussian noise to obtain the data set $`\{x_i\}`$ as follows,

$$x_i=f(\tau _i^{})+\eta \mathrm{𝚛𝚊𝚗𝚍𝚗}(m,1),$$

where $`m=2,3`$ is the dimension of the input space, and randn is Matlab’s standard normal distribution. In Figure 1, in the first row from left to right, we plot the color-coded sample data points corresponding to the following three one-variable functions

$$\begin{array}{ccccc}\hfill f(\tau )& =& (10\tau ,\mathrm{\hspace{0.17em}\hspace{0.17em}10}\tau ^3+2\tau ^210\tau )^T,\hfill & \tau [1,\mathrm{\hspace{0.17em}\hspace{0.17em}1}],\hfill & \hfill \eta =0.1,\\ \hfill f(\tau )& =& (\tau \mathrm{cos}(\tau ),\tau \mathrm{sin}(\tau ))^T,\hfill & \tau [0,\mathrm{\hspace{0.17em}\hspace{0.17em}4}\pi ],\hfill & \hfill \eta =0.2,\\ \hfill f(\tau )& =& (3\mathrm{cos}(\tau ),\mathrm{\hspace{0.17em}\hspace{0.17em}3}\mathrm{sin}(\tau ),\mathrm{\hspace{0.17em}\hspace{0.17em}3}\tau )^T,\hfill & \tau [0,\mathrm{\hspace{0.17em}\hspace{0.17em}4}\pi ],\hfill & \hfill \eta =0.2.\end{array}$$

In the second row, we plot $`\tau _i^{}`$ against $`\tau _i`$, where $`\tau _i`$’s are the computed coordinates by LTSA. Ideally, the $`(\tau _i^{},\tau _i)`$ should form a straight line with either a $`\pi /4`$ or $`\pi /4`$ slope.

It is also important to better understand the failure modes of LTSA, and ultimately to identify conditions under which LTSA can truly uncover the hidden nonlinear structures of the data points. As we have shown in the error analysis in section 6, it will be difficult to align the local tangent information $`\mathrm{\Theta }_i`$ if some of the $`P_i`$’s defined in (3.2) are close to be singular. One effect of this is that the computed coordinates $`\tau _i`$ and its neighbors may be compressed together. To clearly demonstrate this phenomenon, we consider the following function,

$$f(\tau )=[\mathrm{cos}^3(\tau ),\mathrm{sin}^3(\tau )]^T,\tau [0,\pi ].$$

The Jacobi matrix (now a single vector since $`d=1`$) given by

$$J_f(\tau )=1.5\mathrm{sin}(2\tau )[\mathrm{cos}(\tau ),\mathrm{sin}(\tau )]^T$$

is equal to zero at $`\tau =\pi /2`$. In that case the $`\theta `$-vector $`\mathrm{\Theta }_i`$ defined in (4.5) will be computed poorly in the presence of noise. Usually the corresponding $`\mathrm{\Theta }_i`$ will be small which also results in small $`\tau _i`$ and the neighbors of $`\tau _i`$ will also be small. In the first column of Fig 2, we plot the computed results for this 1-D curve. We see clearly near the singular point $`\tau =\pi /2`$ the computed $`\tau _i`$’s become very small, all compressed to a small interval around zero. In the second column of Fig 2, we examine another 1D curve defined by

$$f(\tau )=[10\mathrm{cos}(\tau ),\mathrm{sin}(\tau )]^T,\tau [\pi /2,\mathrm{\hspace{0.17em}\hspace{0.17em}3}\pi /2].$$

We notice that similar phenomenon also occurs near the point $`\tau =\pi `$ where the curvature of the curve is large, the computed $`\tau _i`$’s near the corresponding point also become very small, clustering around zero.

Next we look at the issues of the interaction of sampling density and noise levels. If there are large noises around $`f(\tau _i)`$ relative to the sampling density near $`f(\tau _i)`$, the resulting centered local data matrix $`X_i(I\frac{1}{k}ee^T)`$ will not be able to provide a good local tangent space, i.e., $`X_i(I\frac{1}{k}ee^T)`$ will have singular values $`\sigma _d`$ and $`\sigma _{d+1}`$ that are close to each other. This will result in a nearly singular matrix $`P_i=Q_i^TJ_i`$, and when plotting $`\tau _i^{}`$ against $`\tau _i`$, we will see the phenomenon of the computed coordinates $`\tau _i`$ getting compressed, similar to the case when the generating function $`f(\tau )`$ has singular and/or near-singular points. However, in this case, the result can usually be improved by increasing the number of neighbors used for producing the shifted matrix $`X_i(I\frac{1}{k}ee^T)`$. In Fig 3, we plot the computed results for the generating function

$$f(\tau )=3\tau ^3+2\tau ^22\tau ,\tau [1.1,\mathrm{\hspace{0.17em}1}].$$

The data set is generated by adding noise in a relative fashion,

$$x_i=f(\tau _i)(1+\eta ϵ_i)$$

with normally distributed $`ϵ_i`$. The first three columns in Fig 3 correspond to the noise levels $`\eta =0.01`$, $`\eta =0.03`$, and $`\eta =0.05`$, respectively. For the three data sets, We use the same number of neighbors, $`k=10`$. With the increasing noise level $`\eta `$, the computed $`\tau _i`$’s get expressed at points with relatively large noise. The quality of the computed $`\tau _i`$’s can be improved if we increase the number of neighbors as is shown on the column (d) in Fig 3. The improved result is for the same data set in column (c) with $`k=20`$ used.

As we have shown in Fig 3 (column (d)), different neighborhood size $`k`$ will produce different embedding results. In general, $`k`$ should be chosen to match the sampling density, noise level and the curvature at each data points so as to extract an accurate local tangent space. Too few neighbors used may result in a rank-deficient tangent space leading to over-fitting, while too large a neighborhood will introduce too much bias and the computed tangent space will not match the local geometry well. It is therefore worthy of considering variable number of neighbors that are adaptively chosen at each data point. Fortunately, our LTSA algorithm seems to be less sensitive to the choice of $`k`$ than LLE does as will be shown later.

We now apply LTSA to a $`2`$-D manifold embedded in a 100 dimensional space. The data points are generated as follows. First we generate $`N=5000`$ $`3`$D points,

$$x_i=(t_i,s_i,h(t_i,s_i))^T+0.01\eta _i$$

with $`t_i`$ and $`s_i`$ uniformly distributed in the interval $`[1,\mathrm{\hspace{0.17em}\hspace{0.17em}1}]`$, the $`\eta _i`$’s are standard normal. The $`h(t,s)`$ is a peak function defined by

$$h(t,s)=0.3(1t)^2e^{t^2(s+1)^2}(0.2tt^3s^5)e^{t^2s^2}0.1e^{(t+1)^2s^2}.$$

This function is plotted in the left of Figure 4. We generate two kinds of data points $`x_i^Q`$ and $`x_i^H`$ in 100D space,

$$x_i^Q=Qx_i,x_i^H=Hx_i,$$

where $`Q`$ is a random orthogonal matrix resulting in an orthogonal transformation and $`H`$ a matrix with its singular values uniformly distributed in $`(0,\mathrm{\hspace{0.17em}1})`$ resulting in an affine transformation. Figure 4 plots the coordinates for $`x_i^Q`$ (middle) and $`x_i^H`$ (right).

One advantage of LTSA over LLE is that using LTSA we can potentially detect the intrinsic dimension of the underlying manifold by analyzing the local tangent space structure. In particular, we can examine the distribution of the singular values of the local data matrix $`X_i`$ consisting of the data points in the neighborhood of each data point $`x_i`$. If the manifold is of dimension $`d`$, then $`X_i`$ will be close to a rank-$`d`$ matrix. We illustrate this point below. The data points are $`x_i^Q`$ of the $`2`$D peak manifold in the $`100`$D space. For each local data matrix $`X_i`$, let $`\sigma _{j,i}`$ be the $`j`$-the singular value of the centered matrix $`X_i(I\frac{1}{k}ee^T)`$. Define the ratios

$$\rho _i^{(j)}=\frac{\sigma _{j+1,i}}{\sigma _{j,i}}.$$

In Fig 5, we plot the rations $`\rho _i^{(1)}`$ and $`\rho _i^{(2)}`$. It clearly shows the feature space should be $`2`$-dimensional.

Next, we discuss the issue of how to use the global coordinates $`\tau _i`$’s as a means for clustering the data points $`x_i`$’s. The situation is illustrated by Figure 6. The data set consists of three bivariate Gaussians with covariance matrices $`0.2I_2`$ and mean vectors located at $`[1,1],[1,1],[1,0]`$. There are $`100`$ sample points from each Gaussian. The thick curve on the right panel represents the principal curve computed by LTSA and the thin curve by LLE. It is seen that the thick curve goes through each of the Gaussians in turn, and the corresponding global coordinates (plotted in the middle panel) clearly separate the three Gaussians. LLE did not perform as well, mixing two of the Gaussians.

Last we look at the results of applying LTSA algorithm to the face image data set . The data set consists of a sequence $`698`$ $`64`$-by-$`64`$ pixel images of a face rendered under various pose and lighting conditions. Each image is converted to an $`m=4096`$ dimensional image vector. We apply LTSA with $`k=12`$ neighbors and $`d=2`$. The constructed coordinates are plotted in the middle of Figure 7. We also extracted four paths along the boundaries of of the set of the 2D coordinates, and display the corresponding images along each path. It can be seen that the computed 2D coordinates do capture very well the pose and lighting variations in a continuous way.

## 9 Conclusions and Feature Works

In this paper, we proposed a new algorithm (LTSA) for nonlinear manifold learning and nonlinear dimension reduction. The key techniques we used are construction of local tangent spaces to represent local geometry and the global alignment of the local tangent spaces to obtain the global coordinate system for the underlying manifold. We provide some careful error analysis to exhibit the interplay of approximation accuracy, sampling density, noise level and curvature structure of the manifold. In the following, we list several issues that deserve further investigation.

1. To make LTSA (similarly LLE) more robust against noise, we need to resolve the issue where several of the smallest eigenvalues of $`B`$ are about the same magnitude. This can be clearly seen when the manifold itself consists of several disjoint components. If this is the case, one needs to break the matrix $`B`$ into several diagonal blocks, and apply LTSA to each of the block. However, with noise, the situation becomes far more complicated, several eigenvectors corresponding to near-zero eigenvalues can mix together, the information of the global coordinates seems to be contained in the eigen-subspace, but how to unscramble the eigenvectors to extract the global coordinate information needs more careful analysis of the eigenvector matrix of $`B`$ and various models of the noise. Some preliminary results on this problem have been presented in .

2. The selection of the set of points to estimate the local tangent space is very crucial to the success of the algorithm. Ideally, we want this set of points to be close to the tangent space. However, with noise and/or at the points where the curvature of the manifold is large, this is not an easy task. One line of ideas is to do some preprocessing of the data points to construct some restricted local neighborhoods. For example, one can first compute the minimum Euclidean spanning tree for the data set, and restrict the neighbors of each point to those that are linked by the branches of the spanning tree. This idea has been applied in self-organizing map . Another idea is to use iterative-refinement, combining the computed $`\tau _i`$’s with the $`x_i`$’s for neighborhood construction in another round of nonlinear projection. The rationale is that $`\tau _i`$’s as the computed global coordinates of the nonlinear manifold may give a better measure of the local geometry.

3. A discrete version of the manifold learning can be formulated by considering the data points as the vertices of an undirected graph . A quantization of the global coordinates specifies the adjacency relation of those vertices, and manifold learning becomes discovering whether an edge should be created between a pair of vertices or not so that the resulting vertex neighbors resemble those of the manifold. We need to investigate a proper formulation of the problem and the related optimization methods.

4. From a statistical point of view, it is also of great interest to investigate more precise formulation of the error model and the associated consistency issues and convergence rate as the sample size goes to infinity. The learn-ability of the nonlinear manifold also depends on the sampling density of the data points. Some of the results in non-parametric regression and statistical learning theory will be helpful to pursue this line of research.