text:
\begin{document}
\maketitle

\begin{abstract}
We prove that the pure braid groups on closed, orientable surfaces are
bi-orderable, and that the pure braid groups on closed, non-orientable
surfaces have generalized torsion, thus they are not bi-orderable.
\footnote{\hspace{-.65cm} Keywords: Braid - Surface -
Orderable group.}
\footnote{\hspace{-.65cm} {\em Mathematics Subject Classification:}
 Primary: 20F36. Secondary: 57N05.}
\footnote{\hspace{-.65cm} Partially supported by
DGESIC-PB97-0723 and by the european network TMR Sing. Eq. Diff. et Feuill.}
\end{abstract}

\section{Introduction}

 The purpose of this paper is to answer the following question:
Are pure braid groups on closed surfaces bi-orderable? We will prove
that the answer is positive for orientable surfaces, and negative
for the non-orientable ones.

 In this section we give the basic definitions and
classical results. We also explain what is known about orders on braid
groups, and finally we state our results. In Section~2 we prove the
stated result about closed, orientable surfaces. The
non-orientable case is treated in Section~3.

\subsection{Braids on surfaces}

 Let $M$ be a compact surface, not necessarily orientable,
and let ${\cal P}=\{P_1,\ldots,P_n \}$ be a set of $n$ distinct
points in $M$. Define a {\em $n$-braid based at ${\cal P}$} to be a
collection $b = (b_1,\ldots,b_n)$ of disjoint smooth paths in
$M\times [0,1]$, called {\em strings} of $b$, such that the $i$-th string
$b_i$ runs monotonically in $t\in [0,1]$ from the point $(P_i,0)$ to some
point $(P_j,1)$, $\; P_j\in {\cal P}$.

  An {\em isotopy} is defined as a deformation through braids,
which fixes the ends.  Multiplication of braids is defined by
concatenation. The isotopy classes of braids with this multiplication
form the group $B_n(M,{\cal P})$, called {\em braid group with $n$
strings on $M$ based at ${\cal P}$}.
Note that the group $B_n(M,{\cal P})$ does not depend,
up to isomorphism, on the set ${\cal P}$ of points but only on the
cardinality $n=|{\cal P}|$. So we may write $B_n(M)$ in place of
$B_n(M,{\cal P})$. Further details can be found in \cite{birman}.

 A {\em pure braid} is an element $b\in B_n(M)$ such that $b_i$ ends at
$(P_i,1)$ for all $i=1,\ldots,n$. In other words, $b$ induces the trivial
permutation on ${\cal P}$. Pure braids form a normal subgroup of
$B_n(M)$ called  {\em pure braid group with $n$
strings on $M$}, and denoted by $PB_n(M)$.

 Let $D$ be a closed disc. Notice that $B_n=B_n(D)$ is the classical braid
group of Artin \cite{artin}. Recall that $B_n$
is generated by the set $\{\si_1,\ldots,\si_{n-1}\}$, where $\si_i$ is
the braid shown in Figure~\ref{sigmai}.

\vspace{.3cm}
\begin{figure}[ht]
\centerline{\input{sigmai.pstex_t}}
\caption{The braid $\si_i\in B_n(D)$.}\label{sigmai}
\end{figure}

\subsection{Orderable groups}

 A group $G$ is said to be {\em right-orderable} if there exists a strict
total ordering  $<$ on its elements which is invariant under
right-multiplication: $g<h$ implies $gk < hk$ for all $g,h,k\in G$.
If $<$ is also invariant under left-multiplication, the $G$ is said to be
{\em bi-orderable}.

 The following is a well known characterization of right-orderable
and bi-orderable groups.

\begin{proposition}
 A group $G$ is right-orderable if and only if there exists a subset
${\mathfrak p}\subset G$ such that ${\mathfrak p}^2 \subset {\mathfrak p}$ and
$G= {\mathfrak p} \coprod \{ 1\} \coprod {\mathfrak p}^{-1}$. Moreover, $G$ is
bi-orderable if and only if there exists such ${\mathfrak p}$ also satisfying
$g {\mathfrak p} g^{-1}\subset {\mathfrak p}$, for all $g\in G$.
\end{proposition}

\begin{proof}
 If $G$ is right-orderable (bi-orderable), just take ${\mathfrak p}=
\{g\in G; \; 1<g\}$, the set of positive elements. Conversely, if there
exists ${\mathfrak p}$ verifying the required hypothesis, then define the
order $<$ by: $g<h$ if and only if $hg^{-1}\in {\mathfrak p}$.
\end{proof}

\vspace{.3cm}
 Let us state three properties of orderable groups, which can be found
either in \cite{paris} or in \cite{rolfsenzhu}.

\begin{proposition}\label{proext}
 Let $1 \rightarrow A \stackrel{\al}{\rightarrow} B
\stackrel{\be}{\rightarrow} C \rightarrow 1$ be an exact sequence of groups.
If $A$ and $C$ are right-orderable, with sets of positive elements
${\mathfrak p}_A$ and ${\mathfrak p}_C$ respectively, then the set
${\mathfrak p}_B= \al ({\mathfrak p}_A)\cup \be^{-1}({\mathfrak p}_C)$ defines a
right order on $B$. Moreover, if $A$ and $C$ are bi-orderable, the
order on $B$ defined in such a way is a bi-order if and only if
$\al({\mathfrak p}_A)$ is normal in $B$ (that is, if the order on $A$ is
preserved by conjugation in $B$).
\end{proposition}

\begin{definition}
 We say that a group $G$ has {\em generalized torsion} if there exist
$g,h_1,\ldots,h_k\in G$, $g\neq 1$, such that
$(h_1 g h_1^{-1})(h_2 g h_2^{-1})\cdots (h_k g h_k^{-1})=1$.
\end{definition}

\begin{proposition}\label{progt}
A bi-orderable group has no generalized torsion (in particular, it is
torsion-free).
\end{proposition}

\begin{proposition}
Let $G$ be a right-orderable group, and let $R$ be a ring with no zero
divisors. Then the group ring $RG$ has no zero divisors. Moreover, the
only units of $RG$ are the monomials $rg$, with $r$ invertible in $R$.
\end{proposition}

\subsection{Ordering braid groups}

 Artin braid groups $B_n=B_n(D)$ are known to be right-orderable
(\cite{dehornoy}, see also \cite{fenn}). However they are not
bi-orderable if $n\geq 3$,
since they have generalized torsion \cite{neuwirth}. Nevertheless,
$P_n=PB_n(D)$ is bi-orderable, as shown in \cite{rolfsenzhu} and in
\cite{kimrolfsen}.

 In more generality, it is shown in \cite{rourkewiest} (see also
\cite{shortwiest}) that the
mapping class groups of compact surfaces with boundary, fixing a finite
number of points, are right-orderable. Braid groups on compact surfaces
with boundary are subgroups of these mapping class groups, thus they are
also right-orderable.

 On the other hand, it is not known if the braid groups on closed surfaces
are right-orderable or not. If $M$ is a closed surface different from
the sphere and from the projective plane, then $B_n(D)\subset B_n(M)$
(see \cite{parisrolfsen}). Hence, $B_n(M)$ also has generalized torsion,
and in conclusion, it cannot be bi-orderable.

 In this paper, we will study the pure braid groups of closed surfaces.
We will show the following.

\begin{theorem}\label{teoori}
 If $M$ is a closed, orientable surface, then $PB_n(M)$ is bi-orderable.
\end{theorem}

\begin{theorem}\label{teonon}
 If $M$ is a closed, non-orientable surface, then $PB_n(M)$ has generalized
torsion, for $n\geq 2$. Therefore, $PB_n(M)$ is not bi-orderable
for $n\geq 2$.
\end{theorem}

\section{Closed, orientable surfaces}

 In this section we prove Theorem~\ref{teoori}. In Subsection~\ref{free}
we state that free groups and fundamental groups of orientable
surfaces are bi-orderable. In the case of free groups, we will explicitly
define a bi-order. Then we see in Subsection~\ref{final} that
$PB_n(M)$ is an extension of two groups, $K_n$ and $\pi_1(M)^n$,
which are both bi-orderable. Moreover, the hypothesis of
Proposition~\ref{proext} are satisfied, so $PB_n(M)$ turns to be
bi-orderable.

\subsection{Bi-order of free groups and fundamental groups}\label{free}

 We will explicitly define a bi-order on a given free group using the
so-called {\em Magnus expansion} \cite{magnus}. Let $F$ be a free group
with free system of generators $G=\{x_i \}_{i\in I}$ ($I$ not
necessarily finite). Let $\ZZ\lbra X_I \rbra$ be the ring of formal
power series over the non-commutative indeterminates $\{X_i\}_{i\in I}$.
The {\em Magnus expansion} of $F$ is a multiplicative homomorphism
$M:\: F\rightarrow \ZZ\lbra X_I \rbra$, such that
$ M(x_i)= 1+ X_i$ and  $ M(x_i^{-1})=1-X_i+X_i^2-\cdots$,
for all $i\in I$. $M$ is known to be a well defined and
injective homomorphism, whose image is contained in
$\{1+\eta \in \ZZ\lbra X_I \rbra;\; \eta(0)=0 \}$.

 Let us define a total order on $\ZZ\lbra X_I \rbra$. First, we choose
a total order on the set $\{X_i \}_{i\in I}$. Then, we order the monomials
of $\ZZ\lbra X_I \rbra$ as follows:
$$
  m_1 < m_2 \Leftrightarrow \left\{
\begin{array}{l}
 \mbox{deg}(m_1)< \mbox{deg}(m_2)  \\
 \mbox{deg}(m_1)= \mbox{deg}(m_2)  \mbox{ and } m_1 <_{\mbox{lex}} m_2,
\end{array} \right.
$$
where deg means total degree (the sum of all exponents)  and
$<_{\mbox{lex}}$ means smaller in the lexicographical order.
The total order $\prec$ on $\ZZ\lbra X_I \rbra$ is defined as follows:
given $f,g\in \ZZ\lbra X_I \rbra$, we say that $f\prec g$ if and only if
the coefficient of the smallest non-trivial term of $g-f$ is positive.

 Turning back to $F$, we define the {\em Magnus order} on it: given
$a,b\in F$, we say that $a<b$ if and only if $M(a)\prec M(b)$. The
Magnus order is known to be a bi-order on $F$. Moreover,
one has:

\begin{theorem}{\em \cite{kimrolfsen}}\label{teoh1}
The Magnus order on $F$ is preserved
under any $\Phi\in \mbox{\em Aut}(F)$ which induces the identity on
$H_1(F)=F/[F,F]$. $\fbox{}$
\end{theorem}

\vspace{.3cm}
Let $\psi$ be a permutation of the set $\{X_i \}_{i\in I}$, and consider
its extension $\Psi\in \mbox{\em Aut}(F)$. One has:

\begin{theorem}\label{teopsi}
 If $\psi$ preserves the order on $\{X_i \}_{i\in I}$,
then $\Psi$ preserves the Magnus order on $F$.
\end{theorem}

\begin{proof}
 Notice that, under the action of such a $\Psi$, the degree and the
lexicographical order on the monomials are preserved.
Hence, $\Psi$ preserves the order we defined on the monomials, thus the
order $\prec$ on $\ZZ\lbra X_I \rbra$. Therefore, the Magnus order
on $F$ is also preserved.
\end{proof}

\vspace{.3cm}
 We finish this subsection with the following result.

\begin{theorem}{\em \cite{baumslag}}\label{teopi1}
If $M$ is a closed, orientable surface, then $\pi_1(M)$ is a bi-orderable
group. $\fbox{}$
\end{theorem}

\subsection{$PB_n(M)$ is bi-orderable}\label{final}

 Let $M$ be a closed, orientable surface. Given a pure braid
$b=(b_1,\ldots,b_n)\in PB_n(M)$, we can consider, for all $i=1,\ldots,n$,
the projection $\mu_i$ of $b_i$ over $M$. Since $b\in PB_n(M)$, $\mu_i$
is a loop in $M$ for all $i=1,\ldots,,n$, which represents an element
of $\pi_1(M)$. This defines an epimorphism $\theta:\: PB_n(M)
\rightarrow \pi_1(M)^n$, which sends $(b_1,\ldots,b_n)$ to
$(\mu_1,\ldots,\mu_n)$ (see \cite{birman}).

 Define $K_n=\mbox{ker}(\theta)$. One has the exact sequence
$$
   1 \longrightarrow K_n \longrightarrow PB_n(M)
\stackrel{\theta}{\longrightarrow} \pi_1(M)^n \longrightarrow 1.
$$

\begin{figure}[ht]
\centerline{\input{polairais.pstex_t}}
\caption{The polygon representing $M$ and the braids
$a_{i,{2k+1}}$ and $a_{i,2k}$.}\label{polairais}
\end{figure}

Let us prove that $K_n$  is bi-orderable. First, we represent $M$ as
a polygon of $4g$ sides, identified in the way of Figure~\ref{polairais}.
For all $i=1,\ldots,n$ and all $r=1,\ldots,2g$, we define the braid
$a_{i,r}\in PB_n(M)$ as in Figure~\ref{polairais}:
The $i$-th string of $a_{i,r}$ is
$(s_{i,r}(t),t)\in M \times [0,1]$,
where $s_{i,r}$ is a loop in $M$ based at $P_i$ which goes through the wall
$\al_r$; it goes upwards if $r$ is odd and downwards if $r$ is even.
The $j$-th string of $a_{i,r}$ is $(P_j,t)$ (the trivial string) for all
$j\neq i$.

 Let now $\Om=\{\om_1,\ldots,\om_{2g} \}$ be a set of generators of
$\pi_1(M)$, where $g$ is the genus of $M$. Take $\Om$ in such a way
that $\pi_1(M)=\left< \Om; \: \om_1 \cdots \om_{2g} \om_1^{-1}\cdots
\om_{2g}^{-1}=1 \right>$. For all $\ga\in \pi_1(M)$, choose
a unique word $\widetilde{\ga}$ over $\Om \cup \Om^{-1}$ representing
$\ga$. We denote by $\widetilde{\ga}_{(i)}$ the pure braid obtained
from $\widetilde{\ga}$ by replacing $\om_r^{\pm 1}$ with $a_{i,r}^{\pm 1}$.
Now, for all $i,j\in \{1,\ldots,n \}$, $i<j$, we define the braid
$$
  t_{i,j}=t_{j,i}= \si_i \cdots \si_{j-2} \si_{j-1}^2 \si_{j-2}^{-1}
 \cdots  \si_i^{-1} \in PB_n(M).
$$
Finally, for all $i,j\in \{1,\ldots,n \}$, $i\neq j$, and all
$\ga\in \pi_1(M)$, we define $f_{i,j,\ga}= \widetilde{\ga}_{(i)}
 t_{i,j} \widetilde{\ga}_{(i)}^{-1}$.

\begin{theorem}{\em \cite{gmparis}}\label{teodecomp}
 One has $K_n=( F_n \rtimes(F_{n-1}\rtimes(\cdots(F_3\rtimes F_2)\cdots)))$,
where for all $i=1,\ldots,n-1$, $F_{(n+1)-i}$ is the free group
freely generated by ${\cal F}_{i,n}=\{f_{i,j,\ga};\; i< j\leq n, \:
 \ga\in \pi_1(M)\}$. Moreover, for all $m=2,\ldots,n-1$,
$K_m=(F_m \rtimes(\cdots(F_3\rtimes F_2)\cdots ))$ acts trivially on
$H_1(F_{m+1})$. $\fbox{}$
\end{theorem}

\begin{corollary}\label{cor}
 $K_n$ is bi-orderable.
\end{corollary}

\begin{proof}
 We argue by induction on $n$. If $n=2$, then $K_n=F_2$ is a free group
(of infinite rank), so it is bi-orderable. Suppose then that $n>2$,
and that $K_{n-1}$ is bi-orderable. By Theorem~\ref{teodecomp}, we
have an exact sequence
$$
  1 \longrightarrow F_n \longrightarrow K_n \longrightarrow K_{n-1}
\longrightarrow 1,
$$
where $K_n= F_n \rtimes K_{n-1}$. By definition of bi-order,
conjugation by an element of $F_n$ is an automorphism of $F_n$ which
preserves the Magnus order. We also know, by Theorem~\ref{teodecomp},
that conjugation by an element of $K_{n-1}$ is
an automorphism of $F_n$ which is trivial on $H_1(F_n)$.
Hence, by Theorem~\ref{teoh1}, it also preserves the Magnus
order on $F_n$. Therefore, conjugation by an element of $K_n$ preserves
the Magnus order of $F_n$ and thus, by Proposition~\ref{proext}, $K_n$ is
bi-orderable.
\end{proof}

\vspace{.3cm}
 Let us define an explicit bi-order on $K_n$. First, for all
$i=1,\ldots,n-1$, we order ${\cal F}_{i,n}$ as follows:
$$
  f_{i,j,\ga}< f_{i,k,\de} \Leftrightarrow \left\{
\begin{array}{l}
 j<k  \\   j=k \mbox{ and } \ga <_{\pi_1} \de,
\end{array} \right.
$$
where $<_{\pi_1}$ is a fixed bi-order of $\pi_1(M)$. Then, we consider
the Magnus order on each $F_{(n+1)-i}$ corresponding to this order on
${\cal F}_{i,n}$. The bi-order on $K_n$ which yields from
Corollary~\ref{cor} is the following: for $k,k'\in K_n$, write
$k=k_1k_2\cdots k_{n-1}$ and $k'=k'_1k'_2\cdots k'_{n-1}$,
where $k_i,k'_i\in F_{(n+1)-i}$. Then $k<k'$ if and only
if $k_j<k'_j$ for the greatest $j$ such that $k_j\neq k'_j$.

\vspace{.3cm}
\noindent {\sc Proof of Theorem~\ref{teoori}}: The direct product of
bi-orderable groups is clearly bi-orderable, hence, by
Theorem~\ref{teopi1}, $\pi_1(M)^n$ is bi-orderable. So, by
Proposition~\ref{proext}, we only need to show that conjugation by
an element of $PB_n(M)$ is an automorphism of $K_n$ which preserves
the order.

 Conjugation by an element of $K_n$ preserves the order by definition
of bi-order. Hence, it suffices to show the above claim for the conjugation
by pre-images under $\theta$ of the generators of $\pi_1(M)^n$. A set of
such pre-images is $\{a_{i,r}; \; i=1,\ldots,n,\; r=1,\ldots,2g \}$.
Now, in [G-MP, Lemma~3.15] it is shown that the following relations hold in
$H_1(K_n)$:
$$
    a_{i,r}f_{j,k,\ga}a_{i,r}^{-1}\equiv \left\{
\begin{array}{ll}
 f_{j,k,\ga}  & \mbox{if } i\neq j,k \\
 f_{j,k,(\om_r\ga)} & \mbox{if }  i=j \\
f_{j,k,(\ga \om_r^{-1})} \hspace{1cm} & \mbox{if }  i=k.
\end{array} \right.
$$
We claim that the action of $a_{i,r}$ preserves the Magnus order
on each $F_m$, $m=2,\ldots,n-1$, and hence, it preserves the order on
$K_n$. Clearly, the action of $a_{i,r}$ on $K_n$ is the composition of
an automorphism $\Psi_{i,r}$ which permutes the generators of each $F_m$,
with an automorphism $\Phi_{i,r}$ which is trivial on $H_1(K_n)$.
Therefore, by Theorems \ref{teoh1} and \ref{teopsi}, it suffices
to prove that the permutation induced by $\Psi_{i,r}$ on
${\cal F}_{j,n}$ $(j=1,\ldots,n-1)$ preserves the defined order
on ${\cal F}_{j,n}$.

 Let then $f_{j,k,\ga},f_{j,l,\de}\in {\cal F}_{j,n}$, where
$f_{j,k,\ga} < f_{j,l,\de}$.

\noindent {\bf Case 1:} If $k<l$, then
$\Psi_{i,r}(f_{j,k,\ga})=f_{j,k,\ga'}<f_{j,l,\de'}=\Psi_{i,r}(f_{j,l,\de})$,
where $\ga'$ and $\de'$ are determined by the above relations.

\noindent {\bf Case 2:} If $k=l$ and $\ga <_{\pi_1} \de$, then there are
three possibilities. First, if $i\neq j,k$, one has
$\Psi_{i,r}(f_{j,k,\ga})=f_{j,k,\ga}<f_{j,k,\de}=\Psi_{i,r}(f_{j,k,\de})$.
If $i=j$, one has $\Psi_{i,r}(f_{j,k,\ga})=f_{j,k,(\om_r\ga)}<
f_{j,k,(\om_r\de)} =\Psi_{i,r}(f_{j,k,\de})$, since
$\om_r\ga <_{\pi_1} \om_r\de$ ($<_{\pi_1}$ is a left-order).
Finally, if $i=k$, then
$\Psi_{i,r}(f_{j,k,\ga})=f_{j,k,(\ga\om_r^{-1})}<f_{j,k,(\de\om_r^{-1})}=
\Psi_{i,r}(f_{j,k,\de})$, since $\ga\om_r^{-1} <_{\pi_1} \de\om_r^{-1}$
($<_{\pi_1}$ is a right-order).

 Therefore, $PB_n(M)$ is a bi-orderable group. $\fbox{}$

\section{Closed, non-orientable surfaces}

 We turn now to prove Theorem~\ref{teonon}. Let $M$ be a closed,
non-orientable surface, and let ${\cal P}=\{P_1,\ldots, P_n \}\subset M$.
Then, there exists a submanifold $N$, $\: {\cal P}\subset N \subset M$,
such that $N$ is homeomorphic to a M\"obius strip.

 Consider the subset $C= \RR \times [0,1] \times [0,1] \subset \RR^3$,
and identify $(x,y,t)\sim (x+1,1-y,t)$ for all $x\in \RR$, and all
$y,t\in [0,1]$. One has $(C/\sim) \simeq N\times [0,1]$.
We choose $P_i=(1/2,p_i)\in N$, where $p_i=\frac{i}{n+1}$ for all
$i=1,\ldots,n$.

 Denote by $\Ga=(\ga_1,\ldots,\ga_n)$ the braid on $N$
defined as follows: $\ga_i(t)=(1/2-t,p_i,t)$, for all $i=1,\ldots,n$.
It is the braid represented in Figure~\ref{gade}, for $n=4$.
Notice that $\Ga$ is not a pure braid (going through the wall reverses
the orientation). Denote by $\De$ the following braid on $N$:
$\De=(\si_1\cdots \si_{n-1})(\si_1\cdots \si_{n-2})\cdots(\si_1\si_2)\si_1$.
It is also drawn in Figure~\ref{gade}, for $n=4$.
Remark that we represent $N\times\{ 0\}$ above
$N\times\{ 1\}$ to agree with the usual orientation of braids
(pointing downwards).

\vspace{.3cm}
\begin{figure}[ht]\label{gade}
\centerline{\input{gade.pstex_t}}
\caption{The braids $\Ga, \De \in B_n(N)$.}
\end{figure}

 Now, given a point $p=(x,y,t)\in C/\sim$, we denote by $\overline{p}$
the image of $p$ under the symmetry of plane $y=1/2$. That is,
$\overline{p}=(x,1-y,t)$. In the same way, given a braid $b\in B_n(N)$,
we denote by $\overline{b}$ the image of $b$ under the same symmetry
($\overline{b}_i(t)=\overline{b_i(t)}$).

\begin{lemma}
 For all $b\in B_n(N)$, one has $\Ga b \Ga^{-1}= \overline{b}$.
\end{lemma}

\begin{proof}
 Let $b=(b_1,\ldots,b_n)\in B_n(N)$. For all $i=1,\ldots,n$, we write
$b_i(t)=(\be_i(t),t)$, where $\be_i$ is a path on $N$. We denote by
$\ep$ the permutation induced by $b$ on ${\cal P}$. Consider the
braid  $c=(c_1,\ldots,c_n)=\Ga b \Ga^{-1}$. For all $i=1,\ldots,n$, one has
$$
  c_i(t)= \left\{
\begin{array}{ll}
  (1/2-3t, p_i, t)  &  \mbox{if } t\in [0, 1/3],  \\
  (\be_{n+1-i}(3t-1),t)  &  \mbox{if } t\in [1/3, 2/3],  \\
  (1/2+(3t-2), p_{\ep(n+1-i)}, t) \hspace{1cm} &  \mbox{if } t\in [2/3, 1].
\end{array} \right.
$$
Now consider the isotopy
$H:\: (N\times[0,1])\times [0,1] \rightarrow N\times[0,1]$ given by
$$
    H(x,y,t,s)= \left\{
\begin{array}{ll}
  (x+3ts, y, t)  &  \mbox{if } t\in [0, 1/3],  \\
  (x+s, y, t)  &  \mbox{if } t\in [1/3, 2/3],  \\
  (x+3(1-t)s, y, t) \hspace{1cm} &  \mbox{if } t\in [2/3, 1],
\end{array} \right.
$$
Then $H(c,0)=c$ and $H(c,1)=1\: \overline{b} \: 1 \simeq \overline{b}$.
\end{proof}

\vspace{.3cm}
\noindent {\sc Proof of Theorem~\ref{teonon}}: Take a closed disc $D$,
$\; {\cal P}\subset D\subset N\subset M$. It is well known that, in
$B_n(D)$, $\; \De \si_i \De^{-1}= \si_{n-i}$ for all $i=1,\ldots,n-1$.
So, the same relation holds in $B_n(N)$ and in $B_n(M)$ (every isotopy can
be extended by the identity outside $D$). Hence, for $i=1,\ldots,n-1$,
one has $(\Ga \De)\si_i (\Ga \De)^{-1}=\Ga \: \si_{n-i}\: \Ga^{-1}=
\overline{\si_{n-i}}=\si_i^{-1}$ in $B_n(N)$, thus in $B_n(M)$.

 Therefore, one has $\si_i \left[(\Ga \De)\si_i (\Ga \De)^{-1} \right]= 1$,
and so, $\si_i^2 \left[(\Ga \De)\si_i^2 (\Ga \De)^{-1} \right]= 1$
in $B_n(M)$. It suffices to notice that $\si_i^2$ and $\Ga\De$ are pure
braids, and that $\si_i^2\neq 1$ in $PB_n(M)$ for all $i=1,\ldots,n$,
where $n\geq 2$, to conclude that $PB_n(M)$ has generalized torsion.
$\fbox{}$


TITLE: Difficulties in understanding a conditional probability question from a textbook
QUESTION [2 upvotes]: In the book "Complete Probability & Statistics 1 for Cambridge International AS & A Level" I found the following question:

The homework diaries and completed homework of two Students A and B,
are examined.
There is a probability of 0.4 that Student A does not make a note in
her diary homework set. She always does the homework if it is written
in her diary, but never does the homework if it is not written in.
There is a probability of 0.8 that Student B writes the homework given
in her diary. When she does this, she will do the homework 90% of the
time. If she has nothing written in her diary then she checks with a
friend, who knows what the homework is 50% of the time; Student B
always does the homework if she is told what it is by her friend.
A) Find the probability that Stundent B does her homework on a particular
night.
B) Find the probability that both students do her homework on a
particular night.
C) If one piece of homework was given to a student
but wasn’t done, find the probability that it was given to Student A.

For A) I got $0.8\cdot 0.9+0.2\cdot 0.5=0.82$.
For B) I got $0.82\cdot 0.6=0.492.$
For C) I got: The probability that the homework is not done by both or one student is $1-0.492=0.508$. Therefore we find that $\frac{0.4}{0.508}\approx 0.7874$
But consulting this forum I found that the solution is apparently $\frac{0.382}{0.382+0.108}=0.752$
I do not understand why this is correct. For me, 0.752 is the solution to the question:  If homework of both was checked and one was not done, find the probability that it was Student A's.
Can somebody tell me if I am missing something?

REPLY [3 votes]: The homework diaries and completed homework of two Students A and B,
are examined.
There is a probability of 0.4 that Student A does not make a note in
her diary homework set. She always does the homework if it is written
in her diary, but never does the homework if it is not written in.
There is a probability of 0.8 that Student B writes the homework given
in her diary. When she does this, she will do the homework 90% of the
time. If she has nothing written in her diary then she checks with a
friend, who knows what the homework is 50% of the time; Student B
always does the homework if she is told what it is by her friend.
C) If one piece of homework was given to a student but wasn’t done,
find the probability that it was given to Student A.

I got: The probability that the homework is not done by both or one
student is $1-0.492=0.508$. Therefore we find that
$\frac{0.4}{0.508}\approx 0.7874.$

Firstly, you are interpreting "one student" as "at least one student", whereas the author probably means "exactly one student".
Secondly, you aren't taking into account the given condition that one piece of homework wasn't done.
Now, P(A did but not B) $= [0.6 \times (0.8\times0.1+0.2\times0.5)]=0.108,$
while P(B did but not A) $= [(1-0.6) \times (1-(0.8\times0.1+0.2\times0.5))]=0.328.$
The required probability is P(B did but not A) $\div$  [P(A did but not B) + P(B did but not A)], which equals $$\frac{0.328}{0.328+0.108}=0.752.$$

But consulting this forum I found that the solution is apparently
$\frac{0.382}{0.382+0.108}=0.752.$
I do not understand why this is correct. For me, 0.752 is the solution
to the question:  If homework of both was checked and one was not
done, find the probability that it was Student A's.

Exactly! But, to be exact, "exactly one was not done". -)
By the way, your transcribed answer from that forum contains two identical typos.


text:
\section{Analysis of \brassens{}}
\label{sec:analysis}

In this section, we first discuss the assumptions that guarantee \brassens{} to recover $W$ given the $k$-sparse $r$-separable matrix $M$,
and then discuss the computational complexity of \brassens{}.

\subsection{Correctness}
\label{subsec:correct}

In this section, we show that, given a $k$-sparse $r$-separable matrix $M$, the \brassens{} algorithm provably solves SSNMF, that is, it is able to recover the correct set of indices $\mathcal{J}$ such that ${W = M(:,\mathcal{J})}$, under a reasonable assumption.

Clearly, a necessary assumption for \brassens{} to be able to solve SSNMF is that no column of $W$ is a $k$-sparse nonnegative linear combinations of other columns of $W$, otherwise kSSNPA might set that column of $W$ to zero, hence might not be able to extract it.
\begin{assumption}\label{as:notcombiw}
No column of $W$ is a nonnegative linear combination of $k$ other columns of $W$.
\end{assumption}

Interestingly, unlike the standard separable case (that is, $k=r$), and although it is necessary in our approach with \brassens{},
Assumption~\ref{as:notcombiw} is not necessary in general to be able to uniquely recover $W$.
Take for example the situation of \cref{fig:3intvertaligned}, with three aligned points in the interior of a triangle, so that $r=6$, $m=3$ and $k=2$.
The middle point of these three aligned points is a 2-sparse combination of the other two, by construction.
If there are data points on each segment joining these three interior points and the exterior vertices, the only solution to SSNMF with $r=6$ is the one selecting these three aligned points.
However, Assumption~\ref{as:notcombiw} is a reasonable assumption for SSNMF.

\begin{figure}[t]
    \begin{minipage}{0.47\linewidth}
      \centering
      \input{img/3intvertaligned}
      \caption{There are three interior vertices. One of them ($W_5$) is a combination of the others two.}\label{fig:3intvertaligned}
   \end{minipage}\hfill

    \centering
    \begin{minipage}{0.47\linewidth}
      \centering
      \input{img/badex}
    \caption{There are two interior vertices. One of them ($W_4$) is a 2-sparse combination of data points ($M_1$ and $M_2$).}\label{fig:badex}
   \end{minipage}
\end{figure}

Unfortunately, Assumption~\ref{as:notcombiw} is not sufficient for \brassens{} to provably recover $W$. In fact, we need the following stronger assumption.

\begin{assumption}\label{as:notcombim}
No column of $W$ is a nonnegative linear combination of $k$ other columns of $M$.
\end{assumption}
This assumption guarantees that a situation such as the one shown on  Figure~\ref{fig:badex} where one of the columns of $W$ is a 2-sparse combination of two data points is not possible. In fact, in that case, if \brassens{} picks these two data points before the interior vertex $W_4$ in between them, it will not be able to identify $W_4$ as it is set to zero within the projection step of kSSNPA.

Interestingly, in the standard separable case, that is, $k=r$, the two assumptions above coincide; this is the condition under which SNPA is guaranteed to work.
Although Assumption~\ref{as:notcombim} may appear much stronger than Assumption~\ref{as:notcombiw}, they are actually generically equivalent given that the entries of the columns of $H$ are generated randomly (that is, non-zero entries are picked at random and follow some {continuous distribution}).
For instance, for $m=3$ and $k=2$, it means that no vertex is on a segment joining two data points. If the data points are generated randomly on the segments generated by any two columns of $W$, the probability for the segment defined by two such data points to contain a column of $W$ is zero. In fact, segments define a set of measure zero in the unit simplex.

We can now provide a recovery result for \brassens{}.
\begin{theorem} \label{th:recovery}

Let $M = W H$ with $W = M(:,\mathcal{J})$ be a $k$-sparse $r$-separable matrix so that $|\mathcal{J}| = r$; see \cref{def:sparsesep}.

We have that
\begin{itemize}
    \item If $W$ satisfies \cref{as:notcombim},
    then the factor $W$ with $r$ columns in SSNMF is unique (up to permutation and scaling)  and
    \brassens{} recovers it.

    \item If $W$ satisfies \cref{as:notcombiw}, the entries of $H$ are generated at random (more precisely, the position of the non-zero entries are picked at random, while their values follows a continuous ditribution) and $k < \rank(M)$, then,  with probability one,  the factor $W$ with $r$ columns in SSNMF is unique (up to permutation and scaling)  and
    \brassens{} recovers it.
\end{itemize}
\end{theorem}
\begin{proof}
Uniqueness of $W$ in SSNMF under \cref{as:notcombim} is straightforward: since the columns of $W$ are not $k$-sparse combinations of other columns of $M$, they have to be selected in the index set $\mathcal{J}$. Otherwise, since columns of $W$ are among the columns of $M$, it would not possible to reconstruct $M$ exactly using $k$-sparse combinations of $M(:,\mathcal{J})$. Then, since all other columns are $k$-sparse combinations of the $r$ columns of $W$ (by assumption), no other columns needs to be added to $\mathcal{J}$ which satisfies $|\mathcal{J}|=r$.

Let us show that,
under \cref{as:notcombim}, \brassens{} recovers the correct set of indices $\mathcal{J}$.
kSSNPA can only stop when all columns of $W$ have been identified.
In fact, kSSNPA stops when the reconstruction error is zero,
while, under \cref{as:notcombim}, this is possible only when all columns of $W$ are selected (for the same reason as above).
Then, the postprocessing will be able to identify, among all selected columns, the columns of $W$, because they will be the only ones that are not $k$-sparse combinations of other selected columns.

The second part of the proof follows from standard probabilistic results: since $k < \rank(M)$, the combination of $k$ data points generates a subspace of dimension smaller than that of $\col(M)$. Hence, generating data points at random is equivalent to generating such subspaces at random. Since these subspaces form a space of measure zero in $\col(M)$, the probability for these subspaces to contain a column of $W$ is zero, which implies that \cref{as:notcombim} is satisfied with probability one.
\end{proof}

\subsection{Computational Cost}

Let us derive an upper bound on the computational cost of \brassens{}.
First, recall that solving an NNLS problem up to any precision can be done in polynomial time. For simplicity and because we focus on the non-polynomial part of \brassens{}, we denote $\bar{\mathcal{O}}(1)$ the complexity of solving an NNLS problem.
In the worst case, kSSNPA will extract all columns of $M$. In each of the $n$ iterations of kSSNPA, the problem~\eqref{eq:sparseproj} needs to be solved. When $|\mathcal{J}|=\mathcal{O}(n)$, this requires to solve $n$ times (one for each column of $M$) a $k$-sparse least squares problem in $|\mathcal{J}|=\mathcal{O}(n)$ variables. The latter requires in the worst case $\mathcal{O}(n^k)$ operations by trying all possible index sets; see the discussion after~\eqref{eq:sparseproj}. In total, kSSNPA will therefore run in the worst case in time $\bar{\mathcal{O}}(n^{k+2})$.

Therefore, when $k$ is fixed (meaning that $k$ is considered as a fixed constant) and under \cref{as:notcombim}, \brassens{} can solve SSNMF in polynomial time.
Note that this is not in contradiction with our NP-completeness results when $k$ is fixed (Theorem~\ref{th:nphardssnmf}) because our NP-completeness proof does not rely on \cref{as:notcombim}.

In summary, to make SSNMF hard, we need either $k$ to be part of the input, or the columns of $W$ to be themselves $k$-sparse combinations of other columns of $W$.

config: arxiv

file: 2006.07553/arXiv SSNMF/50analysis.tex


TITLE: Find $x$ if $\tan^{-1}(\frac{1+x}{1-x})=\frac{\pi}{4}+\tan^{-1}x$
QUESTION [1 upvotes]: The relation $\tan^{-1}(\frac{1+x}{1-x})=\frac{\pi}{4}+\tan^{-1}x$ holds true for all $1.$ $ x\in \mathbb R$, $2.$ $ x\in (-1,\infty) $, $3.$ $ x\in (-\infty,1) $, $4.$ $ x\in (-\infty,2)$

I took RHS=$\frac{\pi}{4}+\tan^{-1}x=\tan^{-1}1+\tan^{-1}x=\tan^{-1}(\frac{1+x}{1-x})$,if $x\gt0, x\lt1$
(using the property $\tan^{-1}x+\tan^{-1}y=\tan^{-1}(\frac{x+y}{1-xy})$, if $x\gt0,y\gt0,xy\lt1$)
Therefore, my answer is $x\in(0,1)$, which is not an option.
The answer in the key is $2. $ $x\in(-\infty,1)$.
Looks like, they just used $x\lt1$ condition and not $x\gt0$. Why?

REPLY [0 votes]: Taking the tangent of the two members, you check that this equation holds to a multiple of $\pi$.
The function $\arctan$ maps the real axis to the range $\left(-\dfrac\pi4,\dfrac\pi4\right)$ and is odd. Then $\dfrac\pi4+\arctan x\in\left(0,\dfrac\pi2\right)$ and this excludes the case of $x\ge1$.


TITLE: On characterization of weak convergence.
QUESTION [0 upvotes]: I am reading through Theorem 3.2.5 in Durrett, Probability theory and examples, 4th edition.
How can I prove the following?
Let $ Y_n \rightarrow Y$ a.s.. Then, for any open set $G$, $$\liminf_{n \to \infty} 1_G (Y_n) \geq 1_G (Y)$$ holds.
I know it is equivalent to $$ \liminf_{n \to \infty}1_ {Y_{n}^{-1} (G)}  \geq 1_ {Y^{-1} (G)}  $$ where $Y_{n}^{-1} (G)$, $Y^{-1} (G)$are measurable sets.
Any help is appreciated.

REPLY [1 votes]: Let $\omega$ be any point such that $Y_n(\omega) \to Y(\omega)$. If $\omega \notin Y^{-1}(G)$ the the inequality trivially true since the right side is $0$. If $\omega \in Y^{-1}(G)$ then $Y(\omega) \in G$ and hence $Y_n(\omega) \in G$ for $n$ sufficiently large. Hence  $1_{Y{_n}^{-1}(G)}=1$ for $n$ sufficiently large.. Hence $\liminf 1_{Y{_n}^{-1}(G)} \geq 1=1_{Y^{-1}(G)}$.


TITLE: Order and Degree of Differential Equation
QUESTION [3 upvotes]: I am unable to find the order and degree of the following differential equation:
$x^2(dx)^2 + 2xy dx dy + y^2(dy)^2 - z^2(dz)^2 = 0$
My approach:
Take the term $z^2(dz)^2$ to right hand side and then divide throughout by $(dz)^2$
we get  $$\mathrm{\dfrac{x^2(dx)^2}{(dz)^2} + \dfrac{2xy\ dx \ dy}{(dz)^2} +\dfrac{y^2(dy)^2}{(dz)^2} = z^2}$$
By seeing this, we can deduce that order is 1 and degree 2. But when I am looking at the middle term it contains dx dy/(dz)^2 term which is a term of order 2 and hence degree 1. Well am I right here? Is there any other way of finding the degree and order of this total differential equation? Do suggest, if any.

REPLY [0 votes]: Here order is 1 and degree is 2.
As we can write the above differential equation as
$\mathrm x^2{\left(\dfrac{dx}{dz}\right)^2 + 2xy\left(\dfrac{dx}{dz}\right)\left(\dfrac{dy}{dz}\right) +y^2\left(\dfrac{dy}{dz}\right)^2 = z^2}$


TITLE: Integrating $\ln x$ by parts
QUESTION [6 upvotes]: I am asked to integrate by parts $\int \ln(x) dx$. But I'm at a loss isn't there supposed to be two functions in the integral for you to be able to integrate by parts?

REPLY [3 votes]: Just a good point: Wherever you have $$\int p(x)\ln(x)dx$$ in which $p(x)$ is an integrable function; you can use the integration by parts as follows: $$u=\ln(x),~~dv=p(x)dx$$ such that $$\int udv=uv-\int vdu$$


TITLE: Positive integer powers of real valued functions on real line that satisfies intermediate value property
QUESTION [1 upvotes]: Let $f : \mathbb R \to \mathbb R$ be a function. $f$ is said to have intermediate value property if $f$ takes connected subsets of real line to connected subsets of real line i.e. if $a,x,y \in \mathbb R$ are such that $f(x) < a <f(y)$, then $f(z)=a$ for some $z$ in between $x$ and $y$.
My question is : If $f: \mathbb R \to \mathbb R$ satisfies intermediate value property , then for every integer $n \ge 2$ , does $ f_n :\mathbb R \to \mathbb R$ defined as $ f_n (x)=(f(x))^n,\forall x \in \mathbb R$ , also satisfies intermediate value property ?

REPLY [1 votes]: Notice that if $f:\mathbb{R}\to \mathbb{R}$ and $g:\mathbb{R}\to \mathbb{R}$ satisfy the intermediate value property, then so does $g\circ f$.
Now take $g(x)=x^n$.


TITLE: Find average percent needed over tests to achieve a grade
QUESTION [0 upvotes]: Simple question my mind is blanking on, but I'm trying to work out a formula for finding the average percentage score required across a number of tests given the current average.

So say there are $x$ number of tests remaining,
The student has already done $y$ number of tests of varied results,
If the student is after a final overall pass grade of 50%, what formula could be used to account for their current average and what the average percent they need on their remaining tests to achieve a pass.

E.g, if the student has 7 assessment items, is after a 50% final grade , they achieved 72% on the first assessment, 44% on the second assessment, and 48% on the third assessment, what average percentage do they need from the next 4 assessments to finish with an average of 50%?
Thanks.

REPLY [0 votes]: If they need an average of $50\%$ from $7$ tests, the total of their percentages needs to be $350\%$. So far they have got $164\%$, so they need another $186\%$ from the remaining $4$ tests, which is an average of $46.5\%$.
The calculation is $$\frac{(x+y)\times\text{pass mark}-\text{total so far}}{x}.$$


TITLE: Find expected value of coloured segment
QUESTION [2 upvotes]: Suppose picking $n$ random points from $[0,1]$ interval. For each point we colour segment of least length to it's neighboring point. If length to both neighboring points is the same then we colour to the left neighbour. Find expeceted value of length coloured sements.

REPLY [1 votes]: Add an $n+1$-th point that doesn't affect the colouring. The expected value of the coloured length is the probability that this extra point lies in a coloured segment.
By symmetry, the extra point is the left-most or right-most with probability $\frac2{n+1}$, in which case it's not coloured, and it's the second from the left or right with probability $\frac2{n+1}$, and then it's coloured.
That leaves the cases in which the extra point has second-degree neighbours on either side. Conditional on their position, the extra point and its immediate neighbours are uniformly distributed between them. The probability that either of the extra point's immediate neighbours is the nearest neighbour of the other (remember that the extra point doesn't affect the colouring) is
\begin{align}
&3!\left(\int_0^\frac13\mathrm dx\int_x^\frac{x+1}2\mathrm dy(y-x)+\int_\frac13^\frac12\mathrm dx\int_x^{2x}\mathrm dy(y-x)+\int_\frac12^1\mathrm dx\int_x^1\mathrm dy(y-x)\right)\\
={}&6\left(\frac{19}{649}+\frac{19}{1296}+\frac1{48}\right)\\
={}&\frac7{18}\;.
\end{align}
Thus the expected value of the coloured length is
$$
\frac2{n+1}\cdot0+\frac2{n+1}\cdot1+\left(1-\frac4{n+1}\right)\cdot\frac7{18}=\frac7{18}+\frac49\frac1{n+1}\;.
$$
This is valid for $n\ge3$; for lower $n$ the same approach yields the values that Hagen gave in a comment, $0$ for $n=0$ and $n=1$ and $\frac13$ for $n=2$.
Here's code that checks this result numerically.


text:
\begin{document}

\title{An alternating minimization algorithm for Factor Analysis}

 \author{Valentina Ciccone, Augusto Ferrante and Mattia Zorzi \footnote{Department of Information Engineering, University of Padova}}

\maketitle

\begin{abstract}
The problem of decomposing a given covariance matrix as the sum of a positive semi-definite matrix of given rank and a positive semi-definite diagonal matrix, is considered.
We present a projection-type algorithm to address this problem. This algorithm appears to perform extremely well and is extremely fast even when the given covariance matrix has a very large dimension. The effectiveness of the algorithm is assessed through simulation studies and by applications to three real datasets that are considered as benchmark for the problem. A local convergence analysis of the algorithm is also presented.
\end{abstract}

\section{Introduction}\label{Intro}

The problem of decomposing a given covariance matrix into the sum of a low rank matrix $L$ plus a diagonal matrix $D$ bursts more than a century of tradition in scientific literature. In fact, it may be viewed as a linear algebraic counterpart of a \textit{Factor Analysis} problem which is a problem in multivariate statistics {\vale aiming to extract statistical commonalities among data.\\ Factor models were first introduced by Spearman \cite{spearman_1904} at the beginning of the last century in the framework of psychological tests. Since then a rich stream of literature has followed combining psychology and mathematics, see for example \cite{BURT_1909}, \cite{Spearman-Holzinger-24}, \cite{kaiser1958varimax}, \cite{tucker1973reliability}, \cite{shapiro1982rank} and \cite{Bekker-deLeeuw_1987} with reference therein. Soon the interest for this type of models has grown significantly also outside the psychological community, see \cite{deistler2007}, \cite{scherrer1998structure},\cite{GEWEKE_DYNAMIC}, \cite{DEISTLER_1997}, \cite{picci1986dynamic}, \cite{bottegal2015modeling} and \cite{zorzi2016ar}, boasting nowadays applications in countless fields of science.
We refer to \cite{finesso2016factor}, \cite{ning2015linear} and \cite{bertsimas2017certifiably} }
for an {\vale up to date} discussion on the importance of the problem, on its applications, on the formidable stream of literature produced on this topic in the last century, and on the numerous variants in which the problem can be formulated.

In this work we take an optimization-oriented viewpoint: for a given covariance matrix $\Sigma$ and a given rank $r$ we want to find a positive semidefinite matrix $L$ with rank not larger than  $r$ and a positive semidefinite diagonal matrix $D$ such that their sum is as closest as possible to $\Sigma$.
A closed-form solution for this problem appears to be out of reach so that numerical techniques are needed. We propose an easy-to-implement iterative algorithm, based on alternating minimization, to solve numerically the considered problem.
This algorithm appears to perform extremely well and in simulations converges very rapidly to the solution.
However, despite the simplicity of the algorithm, the  convergence analysis is non trivial due to the non-convexity of the set of low rank matrices.

The paper is organized as follows. In Section \ref{Section_FA} we first introduce the Factor Analysis problem together with a brief review of the available literature. Then  the  addressed problem is stated and motivated.
In Section \ref{Section_algorithm} we present the proposed algorithm and  the results of numerical simulations are summarized in Section \ref{sim} together with an application to real datasets.  In Section \ref{Section_convergence} the local convergence analysis for the proposed algorithm is discussed.
Finally, Section \ref{pojection-approach} proposes  a different interpretation of the proposed algorithm that can be viewed as an alternating projection procedure.

\subsection{Notation}
Given a matrix $M$, we denote its transpose by $M\tp$; if $M$ is a square matrix
$\tr(M)$ denotes its trace.

The symbol $\mathbf{Q}_n$ denotes the vector space of real symmetric matrices of size $n$. If $X\in\mathbf{Q}_n$ is positive definite or positive semi-definite we write $X\succ 0$ or $X\succeq 0$, respectively.
Moreover, we denote by $\mathbf{D}_n$ the vector space of  diagonal matrices  of size $n$.

We denote by $\ofd(\cdot)$ the self-adjoint operator orthogonally projecting $\mathbf{Q}_n$ onto  the orthogonal complement of $\mathbf{D}_n$ in $\mathbf{Q}_n$, i.e. if $M\in  \mathbf{Q}_n$, $\ofd(M)$ is the matrix in which each off-diagonal element is equal to the corresponding element of $M$ (and each diagonal element is clearly zero).
The Frobenius norm is denoted by $\Vert \cdot \Vert_F$ while $\Vert \cdot \Vert$ denotes the Euclidean norm.

\section{Preliminaries in Factor Analysis and Problem Definition}\label{Section_FA}

Factor models are used to described high dimensional vectors of data in terms of a small number of common latent factors.
In its simplest formulation, the classic (linear static) factor model is given by
 \alg{\label{fact_model} y&=Ax+z}
 where $A\in \mathbb{R}^{n\times r}$, with $r<<n$, is the so-called factor loading matrix, $x$ is the vector of (independent) latent factors and $z$ represents the idiosyncratic component. Here, $x$ and $z$  are zero-mean, independent Gaussian random vectors; the covariance matrix of $x$  is the identity matrix of dimension $r$ and the covariance matrix of $z$ is a diagonal matrix $D\in \mathbf{D}_n$. Note that, $Ax$ represents the latent variable.
Clearly, $y$ is itself a Gaussian random vector with zero mean and we denote by $\Sigma$ its covariance matrix. Since $x$ and $z$ are independent it holds that
\begin{equation}
\label{decomposition}
\Sigma=L+D
\end{equation}
where $L:=AA^\top$ and $D$ are the covariance matrices of $Ax$ and $z$, respectively. Thus, $L$  has rank equal to $r$, and $D$ is diagonal.

Hence, in its original conception the construction of a factor model is mathematically equivalent to a matrix additive decomposition problem which seeks, for a given $\Sigma$, a decomposition of the type of \eqref{decomposition}.
Of course the model is maximally parsimonious if the rank of $L$ is minimum.
The problem of minimizing the rank of $L$ in  decomposition (\ref{decomposition})  is known as
Frisch's problem and, to date, no exact solution for such a problem is actually available, with the only exception of the special case when this minimum rank is $r=n-1$, in which case a closed-form parametrization of the solutions is provided in \cite{reiersol1950identifiability}. This lack of explicit formulas has motivated a rich stream of literature and different numerical approaches which have been proposed over the years.
A relaxation of this problem has also been considered in which the matrix $D$ is only required to be diagonal but not positive semi-definite. This is known as Shapiro's problem.

The main difficulty in these problems is related to the non-convexity of the rank function so that a viable alternative is to consider the so called  \textit{minimum trace factor analysis} problem, \cite{shapiro1982rank}, \cite{MIN_RANK_SHAPIRO_1982}:
\begin{equation}
\label{classic}
\begin{aligned}
\min_{L,D\in\mathbf{Q}_n} \;\; & \tr (L)\\
& L,D\succeq 0 \\
& \Sigma = L+D\\
& D\in \mathbf{D}_n
\end{aligned}
\end{equation}
 where the trace of L is used as convex surrogate of the rank function as shown in \cite{FAZEL_MIN_RANK_APPLICATIONS_2002}, \cite{FAZEL_MINIMUM_RANK_2004}.

Note that, in many cases the equality constraint in \eqref{classic} may be too compelling. Therefore, an alternative approach is to allow for residuals in the decomposition. Typically, this leads to an optimization problem where the residual $\Sigma - L - D$ is minimized with respect to a chosen norm under a constraint limiting the rank of $L$. This approach is known as \textit{minimum residual factor analysis}, see \cite{harman1966factor}, \cite{shapiro2002statistical},\cite{bertsimas2017certifiably}.
Note that the presence of the rank constraint makes such problems non convex and several heuristic have been proposed to deal with it.\\

Other approaches to factor analysis encompass: principal component factor analysis as in \cite{bai2008large}, maximum likelihood methods as in \cite{bai2012statistical}, or the establishing of a certificate of optimal low rank as in \cite{guttman1954some}.
Moreover, several variants of the mentioned approaches have been proposed by weakening modelling assumptions or by introducing additional constraints for example to account for errors in the covariance matrix estimation as in  \cite{ciccone2017factor}, \cite{ning2015linear} and \cite{agarwal2012noisy}.

The problem we are going to consider is a minimum residual type problem:
for a given $r$ and a given matrix $\Sigma$ we want to find a positive semidefinite matrix $L$ with rank at most $r$ and a positive semidefinite diagonal matrix $D$ such that their sum is as close as possible to $\Sigma$.
This can be formalized as follows:
\begin{equation}
\begin{aligned}
\label{the_problem0}
(L^*,D^*):=\argmin_{L\in \mathcal{L}_{n,r},D\in \mathcal{D}_n} \quad & \Vert \Sigma -L-D \Vert_F^2
\end{aligned}
\end{equation}
where the sets $\mathcal{D}_n$ and $\mathcal{L}_{n,r}$ are defined as:
\begin{equation*}
\mathcal{L}_{n,r}:=\lbrace X\in\mathbf{Q}_n: X\succeq 0, \; \text{rank}(X) \leq r \rbrace,
\end{equation*}
\begin{equation*}
\mathcal{D}_n:=\lbrace X \in\mathbf{D}_n: X\succeq 0\rbrace.
\end{equation*}
 Note that, in practice, $r$ can be obtained by resorting to available methods for estimating the number of factors as \cite{bai2002determining}, \cite{lam2012factor} and \cite{ciccone2017factor}. Alternatively, the problem can be solved for increasing values of $r$ until the residue $\Vert \Sigma -L^*-D^* \Vert_F$ is not greater than a certain tolerance.  In the case that $\Sigma$ is the sample covariance matrix estimated from the data, this is equivalent to find a good trade-off between the fit term (i.e. the residue) and the complexity of the model (i.e. $r$).

{\vale
Our approach is close in spirit to the one proposed by } \cite{bertsimas2017certifiably} where the q-norm of the residue is minimized under the following ulterior constraint:
\beq\label{additionalconstr}
\Sigma-D \succeq 0.
\eeq
This constraint is perfectly justified if we assume that the covariance matrix  $\Sigma$ of
$y$ is known with great precision, $r$ is the number of the most significant common factors and the residue $\Sigma -L^*-D^* $ accounts for other common factors that are less significant.
Our approach considers instead the case (that is realistic in many practical situations) in which $\Sigma$ has been estimated form the data and is therefore only an approximation of the ``true'' covariance matrix; for more details on this case we refer the reader to  \cite{bai2002determining}, \cite{lam2012factor}, \cite{ciccone2017factor} and \cite{fan2013large}. In this setting, the residue $\Sigma -L^*-D^* $ accounts also for
the uncertainty in the estimation of $\Sigma$ so that the constraint $\Sigma-D \succeq 0$ must  not
be imposed.
Of course, if we find an exact decomposition $\Sigma=L^*+D^* $ so that  the residue vanishes, the constraint $\Sigma-D^*  \succeq 0$ is automatically satisfied.

\section{The Proposed Algorithm}\label{Section_algorithm}
A closed-form solution for Problem \eqref{the_problem0} appears to be out of reach.  However, this Problem appears to be well suited for a coordinate descent type iterative algorithm.
Such algorithm alternates between solving a minimization problem with respect to $L$ and  a minimization problem with respect to $D$ in the following fashion:
\begin{equation}
\label{alternate_min}
\begin{aligned}
L_k &= \argmin_{L\in\mathcal{L}_{n,r}} \Vert {\Sigma} - L- D_{k-1}\Vert_F\\
D_k &= \argmin_{D\in\mathcal{D}_n} \Vert {\Sigma} - L_{k}- D\Vert_F
\end{aligned}
\end{equation}
 where $L_k$ and $D_k$ denotes the value of $L$ and $D$, respectively, at the $k$-th iteration. Both these sub-problems admit explicit solutions which are provided by the projection operators onto the set $\mathcal{L}_{n,r}$ and $\mathcal{D}_n$, respectively, as described below.

Let $X\in \mathbf{Q}_n$ and consider its spectral decomposition $X=USU^\top$, $U\in\mathbf{O}_n$ and $S=\diag(s_1,\,...,\, s_n)$ with $s_1\geq s_2 \geq ... \geq s_n$ being the eigenvalues of $X$ arranged in decreasing order. Then, the closest   matrix with rank at most $r$ to $X$ in  the Frobenius norm is obtained applying the projector $P_{\mathcal{L}_{n,r}}$:
\begin{equation*}
P_{\mathcal{L}_{n,r}}(X):=U \diag (f_l (s_1), \, ..., \, f_l(s_n) )U^\top
\end{equation*}
with $f_l(\cdot)$ defined as
\begin{equation}
f_l(s_i):=
\begin{cases}
s_i & \text{for } i\leq r \wedge s_i>0 \\
0 & \text{otherwise.}
\end{cases}
\end{equation}
On the other hand, the projector $P_{\mathcal{D}_n}$ onto the set $\mathcal{D}_n$ is:
\begin{equation}\label{proiezionesud}
P_{\mathcal{D}_n}(X):= \diag(f_d(X_{11}), \, ...., \, f_d(X_{nn}))
\end{equation}
with $f_d(\cdot)$ defined as
\begin{equation}
f_d(X_{ii}):=
\begin{cases}
X_{ii} & \text{if } X_{ii}>0 \\
0 & \text{otherwise.}
\end{cases}
\end{equation}
Then, at $k$-th iteration the algorithm computes:
\[L_k= P_{\mathcal{L}_{n,r}}({\Sigma}-D_{k-1})\]
\[D_k= P_{\mathcal{D}_n}({\Sigma}-L_k).\]
 The complete procedure is outlined in Algorithm 1: $\varepsilon>0$ is the maximum error allowed in the relative decomposition error, while $N$ represents the maximum number of iterations.
\vspace{0.2cm}
\hrule
\vspace{0.1cm}\noindent
\textbf{Algorithm 1 }
\hrule
\vspace{0.2cm}
\textbf{Input:} ${\Sigma}$, $r$, $\epsilon$, N\\
\indent\textbf{Output:}  $L^*, \;D^*$\\
\indent\textbf{Initialize:} initialize $D$ randomly, i=0\\
\indent\textbf{while} $\Vert {\Sigma} -L-D \Vert_F^2/ \Vert {\Sigma} \Vert_F^2 < \epsilon$ \textbf{and} i $<$ N\\
\indent \indent $L= P_{\mathcal{L}_{n,r}}({\Sigma}-D)$\\
\indent \indent $D= P_{\mathcal{D}_n}({\Sigma}-L)$\\
\indent \indent i=i$+1$\\
\indent\textbf{end while} \\
\indent \indent  $L^*=L$, $D^*=D$\\
\indent\vspace{0.05cm}
\hrule
\vspace{0.5cm}

\section{Numerical Simulations}\label{sim}

To provide empirical evidence of the convergence properties of the algorithm simulations studies have been performed
 by using the software Matlab-R2012b on a 2014 laptop MacBook Pro, Quad-i7 2.0 GHz.

 To begin with, we have considered the case of a covariance matrix, $\Sigma$, computed as the sum of a
randomly generated positive semidefinite low-rank matrix $L$ of dimension $n$ and rank $r$, and a randomly generated positive definite diagonal matrix $D$.
We have performed $200$ Monte Carlo runs with $n=40$ and $r=4$ and  $200$  runs with $n=40$ and $r=10$. The original low-rank and diagonal matrices are recovered with negligible numerical errors. Indeed the following quantities:
\begin{itemize}
\item the  relative decomposition error  on $L+D$: $\Vert \Sigma- L^*- D^* \Vert / \Vert \Sigma\Vert$;
\item the relative error on $L$: $\Vert L-L^*\Vert / \Vert L\Vert$;
\item the relative error on $D$: $\Vert D-D^*\Vert / \Vert D\Vert$;
\end{itemize}
\noindent are all of the order of $10^{-10}$. The average computational time for each experiment is less than five hundredths of a second: in less than half minute all $2\times 200$  runs converged.

To account for how the algorithm scales with the dimensionality of the problem two further numerical experiments has been conducted.
First, we have considered the case of a fixed rank, $r=8$, and increasing dimensions: $n=20*2^{j}$, with  $j=0\ldots 5$. For each value of $n$, $50$ factor models have been generated and the resulting covariance matrices serve as input for the algorithm. The statistics of the execution time (in seconds) are summarized in Table \ref{table1}. \\
Second, we have considered the case of a fixed $r/n$ ratio of $0.2$, with $n$ taking values $n=20*2^{j}$, with  $j=0\ldots 5$. For each of them, $50$ factor models have been generated and the resulting covariance matrices serve as input for the algorithm. The statistics of the execution time (in seconds) are summarized in Table \ref{table2}. \\
Both experiments provide evidence that the algorithm scales extremely well with dimensionality.

\begin{table}
\centering

\parbox{.45\linewidth}{
\centering
\begin{tabular}{|l| c| c | c|}
\hline
 $n$ & $r$ & mean & st. dev.\\
\hline
 20 & 8 & 0.0610 & 0.0545\\
 40 & 8 & 0.0349 & 0.0092\\
 80 & 8 & 0.0642 & 0.0081\\
 160 & 8 & 0.1927 & 0.0173\\
 320 & 8 & 1.2286 & 0.0875\\
 640 & 8 & 5.4973 & 0.2793\\
 1280 & 8 & 26.3725 & 0.9813\\
\hline
\end{tabular}
\caption{\small For each value of $n$ the table displays the mean execution time (in seconds) and standard deviation across $50$ experiments.
}
\label{table1}
}
\quad
\parbox{.45\linewidth}{
\centering
\begin{tabular}{|l| c| c | c|}
\hline
 $n$ & $r$ & mean & st. dev.\\
\hline
 20 & 4 & 0.0219 & 0.0452\\
 40 & 8 & 0.0362 & 0.0093\\
 80 & 16 & 0.1069 & 0.0171\\
 160 & 32 & 0.3842 & 0.0444\\
 320 & 64 & 2.7974 & 0.1565\\
 640 & 128 & 14.7181 & 0.6733\\
 1280 & 256 & 85.5031 & 26.8018\\
\hline
\end{tabular}
\caption{\small For each couple of $(n,r)$ the table displays the mean execution time (in seconds) and standard deviation across $50$ experiments.}
\label{table2}
}
\end{table}

 Finally, we have considered the case of a covariance matrix which admits only approximately a "low-rank plus diagonal" decomposition. This case is of practical interest in factor analysis because typically only an estimate, $\hat{\Sigma}_N$, of $\Sigma$ is available.\\
Given a covariance matrix $\Sigma$ generated as before (which therefore admits an exact "low-rank plus diagonal" decomposition), we have generated a sample of numerosity $N$ from the distribution $\mathcal{N}(0, \Sigma)$ and we have estimated the corresponding sample covariance $\hat{\Sigma}_N$ which serves as input for the algorithm. We have considered the same setting as before with $n=40$, $r=4$ and $n=40$, $r=10$. In both cases for each sample size $N= 200, 500, 1000$ we have performed $200$ Monte Carlo runs.
The $6\times 200$ simulations took less than 5 minutes to converge  and we observed the following:
\begin{enumerate}
\item
In all the $6\times 200$ simulations the sequence ($D_k,L_k$) produced by Algorithm 1 converged to a stationary point
($L^*,D^*$)
and, as discussed in Proposition \ref{propconv} below, this point is a (at least) local minimum of the cost function.

\item
In all the $6\times 200$ simulations, the inequality $$\Vert L^*+D^* -\hat{\Sigma}_{N} \Vert_F - \Vert\Sigma_{true} -\hat{\Sigma}_{N} \Vert_F\leq 0$$ is satisfied which provides a sanity check on the performance of the proposed algorithm.
In fact, especially for $N=1000$, $\Sigma_{true}$ may be viewed as a good approximation of $\hat{\Sigma}_{N}$ and, on the other hand, we know that, by construction, $\Sigma_{true}$ may be decomposed as the sum of a low rank positive semidefinite matrix and a diagonal positive matrix. Hence, $\Sigma_{true}=L_{true}+D_{true}$ may be viewed as a benchmark which is always outperformed by the decomposition provided by the proposed algorithm.
\end{enumerate}

The results for the decomposition error $\Vert \hat{\Sigma}_N- L^*- D^* \Vert / \Vert \hat{\Sigma}\Vert$ are summarized in Figures \ref{fig:figure_5} and \ref{fig:figure_6}. Figures \ref{fig:figure_7} and \ref{fig:figure_8} display the following quantities:
\begin{itemize}
\item the relative  decomposition error on $L+D$: $\Vert \Sigma- L^*-D^*\Vert/ \Vert \Sigma \Vert$;
\item the relative error on $L$: $\Vert L-L^*\Vert/ \Vert L\Vert$;
\item the relative error on $D$: $\Vert D-D^*\Vert/ \Vert D\Vert$. \end{itemize}

The obtained results appear extremely promising.
\begin{figure}[h!]
	\centering
	\includegraphics[trim={5cm 2cm 4cm 1cm},clip,width=0.5\textwidth]{03042018_r4_noise_deco_error}
	\caption{ Case $r=4$. Decomposition errors $\Vert \hat{\Sigma}_N- L^*-D^*\Vert/ \Vert\hat{\Sigma}_N \Vert$ with $N=200$, $N=500$ and $N=1000$. }
\label{fig:figure_5}
\end{figure}
\begin{figure}[h!]
	\centering
	\includegraphics[trim={5cm 2cm 4cm 1cm},clip,width=0.5\textwidth]{03042018_r10_noise_deco_error}
	\caption{Case $r=10$. Decomposition errors $\Vert \hat{\Sigma}_N- L^*-D^*\Vert/ \Vert \hat{\Sigma}_N \Vert$  with $N=200$, $N=500$ and $N=1000$.}
\label{fig:figure_6}
\end{figure}
\begin{figure}[h!]
	\centering
	\includegraphics[trim={4.5cm 0.5cm 3.5cm 0.5cm},clip,width=0.7\textwidth ,height=11cm]{03042018_r4_noise_part}
	\caption{ Case $r=4$. The displayed quantities are: $\Vert \Sigma- L^*-D^*\Vert/ \Vert \Sigma \Vert$, $\Vert L-L^*\Vert/ \Vert L\Vert$ and $\Vert D-D^*\Vert/ \Vert D\Vert$,  where $L^{*}$ and $D^{*}$ represent the estimates with $N=200$, $N=500$ and $N=1000$.}
\label{fig:figure_7}
\end{figure}
\begin{figure}[h!]
	\centering
	\includegraphics[trim={4.5cm 0.5cm 3.5cm 0.5cm},clip,width=0.7\textwidth ,height=11cm]{03042018_r10_noise_part}
	\caption{Case $r=10$. The displayed quantities are: $\Vert \Sigma- L^*-D^*\Vert/ \Vert \Sigma \Vert$, $\Vert L-L^*\Vert/ \Vert L\Vert$ and $\Vert D-D^*\Vert/ \Vert D\Vert$, where $L^{*}$ and $D^{*}$ represent the estimates with $N=200$, $N=500$ and $N=1000$. }
\label{fig:figure_8}
\end{figure}

\newpage

\vspace{1cm}

\subsection{Application to real data }
In this section we investigate the performance of the proposed method on three real world datasets which are popular benchmark in factor analysis:
\begin{itemize}
\item the bfi dataset, from the R library psych, which consists of $2800$ observations on 28 variables: 25 variables represent personality self-reported items and while 3 variables represent demographic variables;
\item the neo dataset, also from the R library psych, which consists of a correlation matrix of size $30 \times 30$ estimated from 1000 observations;
\item the Harman dataset, from the R library datasets, which consists of a correlation matrix of size $24 \times 24$ estimated from 145 observations: the cross-section represents psychological tests carried out to seventh- and eighth-grade children.
\end{itemize}
These datasets have been used in \cite[Section 5.3]{bertsimas2017certifiably} to compare the performance of their  approach, which minimizes the $q$-norm of the residue (with $q=1$), against different factor analysis methods. This approach can be considered as the state of the art as it outperforms the other available methods.
In this section we take it as benchmark for comparisons and we repeat the analysis in \cite[Section 5.3]{bertsimas2017certifiably}.\\
The adopted measure of performance is the explained variance, defined as $$\sum_{i=1}^r \lambda_i (L^*)/ \sum_{i=1}^n | \lambda_i (\Sigma - D^*)|.$$
For each dataset Problem \eqref{the_problem0} is solved for the values of $r$ considered in \cite{bertsimas2017certifiably}. The results are depicted in Figure \ref{fig:figure_9}.  The proposed method provides a higher amount of explained variance with respect to the method proposed in \cite{bertsimas2017certifiably} that can be considered to be the state-of-the-art as, so far, it outperforms all the available methods. Moreover, our method  shows a flexibility in delivering different models with varying $r$ which is similar to that of the method proposed in \cite{bertsimas2017certifiably}.

\begin{figure}[h!]
	\centering
	\includegraphics[trim={2cm 0cm 3cm 0cm},clip,width=1\textwidth]{comparison_jml}
	\caption{Proportion of variance explained by the proposed method (blue line) and the benchmark method (red line).}
\label{fig:figure_9}
\end{figure}

In the cases analyzed in these examples the covariance matrix is estimated from a relatively small number of data so that it is reasonable to assume that the residues are not only caused by the presence of less significant latent factors. Thus, we are in the typical situations where our method applies.

\section{Convergence analysis}\label{Section_convergence}
In this section we discuss  the convergence of the proposed algorithm to a local minimum.
First of all we observe that the iterative minimization in \eqref{alternate_min}
produces a sequence of values  for the objective function that is monotonically non-increasing.
Since the objective  function is clearly bounded from below we have the following obvious result.

\begin{lemm}\label{lemseqfob}
For $h\in\mathbb{N}$, define the sequence  $F_h$ by $F_{h}:=\Vert \Sigma -L_{k}-D_k \Vert_F^2 $
for $h=2k$ (even), and $F_{h}:=\Vert \Sigma -L_{k+1}-D_k \Vert_F^2 $, for $h=2k+1$ (odd), where $L_k$, $D_k$ is the sequence produced by Algorithm 1.
Then the sequence $F_h$ is monotonically non-increasing and has limit as $h\rightarrow \infty$.
\end{lemm}

Establishing the convergence for $L_k$ and $D_k$ is less trivial.
We start with $D_k$. To this aim we observe that
as a consequence of Lemma \ref{lemseqfob}, we have that
 $\varepsilon_k:=F_{2k-1}-F_{2k}$   not only converges to zero but it converges sufficiently fast.

\begin{lemm}\label{lemmaalphaauriga}
Assume that  $\varepsilon_k:=F_{2k-1}-F_{2k}$ tends to zero faster than $1/k^{2q}$ with $q>1$ and
let $D_k$ be the sequence of diagonal matrices produced by Algorithm 1.
Then the sequence $D_k$ converges to a certain diagonal matrix $D^*\in \mathcal{D}_n$.
\end{lemm}
\proof
We have
$$F_{2k}=\Vert \Sigma -L_{k}-D_{k} \Vert_F^2=F_{2k-1}-\varepsilon_k=\Vert \Sigma -L_{k}-D_{k-1} \Vert_F^2-\varepsilon_k.$$
Let $s_k(i):=[\Sigma -L_{k}]_{ii}$ be the $i$-th element in the diagonal of $\Sigma -L_{k}$
and $d_k(i):=[D_k]_{ii}$ be the $i$-th element in the diagonal of $D_k$.
Since in (\ref{proiezionesud}) for each $i$, $d_k(i)$ is chosen independently of the others in order to minimize
$\Vert \Sigma -L_{k}-D_{k} \Vert_F^2$, we have that
\bea\nn
-\varepsilon_k&=&\Vert \Sigma -L_{k}-D_{k} \Vert_F^2-\Vert \Sigma -L_{k}-D_{k-1} \Vert_F^2\\
\nn
&=&  \sum_{i=1}^n \lbrace [s_k(i)-d_k(i)]^2-[s_k(i)-d_{k-1}(i)]^2 \rbrace \\
\nn
&\leq &
 [s_k(i)-d_k(i)]^2-[s_k(i)-d_{k-1}(i)]^2
 \eea
 which yields
$$\varepsilon_k \geq [d_{k-1}(i)-d_k(i)][d_k(i)+d_{k-1}(i)-2s_k(i)]. $$
Now, we can consider two cases: if $s_k(i)\geq 0$, then the minimizer $d_k(i)$ is equal to
$s_k(i)$, so that we have
$$\varepsilon_k \geq [d_{k-1}(i)-d_k(i)]^2.$$
If $s_k(i)< 0$, then $d_k(i)=0$  so that we have again
$$\varepsilon_k \geq d_{k-1}(i)[d_{k-1}(i)-2s_k(i)]\geq [d_{k-1}(i)-d_k(i)]^2. $$
In conclusion, in both cases,
we have
$$|d_{k-1}(i)-d_k(i)|\leq \alpha_k:=\sqrt{\varepsilon_k}.$$
As a consequence, we have
\bea\nn
|d_{k+m}(i)-d_k(i)|&\leq & |d_{k+m}(i)-d_{k+m-1}(i)|
+ |d_{k+m-1}(i)-d_{k+m-2}(i)|+\dots \\
\nn
&&
+ |d_{k+1}(i)-d_{k}(i)|\\
\nn
&\leq & \alpha_{k+m}+\dots + \alpha_{k+1}\\
\nn
& = & \sum_{l=1}^m\alpha_{k+l}\\
\nn
&\leq &
\sum_{l=1}^m \frac{M}{(k+l)^{q}}\\
\nn
&\leq &
\sum_{l=1}^\infty \frac{M}{(k+l)^{q}}\\
\nn
&\leq &
 \sum_{h=k+1}^\infty \frac{M}{h^q}
= \left[\sum_{h=1}^\infty \frac{M}{h^q}-\sum_{h=1}^{k+1} \frac{M}{h^q}\right]
\eea
where  $M$ is a constant and $q>1$ so that all the infinite sums converge to a finite value.
Since we have
\begin{align*}
\lim_{k\rightarrow \infty} \left[\sum_{h=1}^\infty \frac{M}{h^q}-\sum_{h=1}^{k+1} \frac{M}{h^q}\right]=\sum_{h=1}^\infty \frac{M}{h^q}-\lim_{k\rightarrow \infty} \sum_{h=1}^{k+1} \frac{M}{h^q}=0,
\end{align*}
we can conclude that $\lim_{l,k\rightarrow\infty}|d_l(i)-d_k(i)|=0$, so that $d_k(i)$ is a Cauchy sequence and hence  it converges.
Since this holds for each $i=1,\dots,n$, we have that the sequence $D_k$ converges to a certain diagonal matrix $D^*$.
Finally, since $\mathcal{D}_n$ is closed, clearly $D^*\in \mathcal{D}_n$.
\qed

For the convergence of the sequence $L_k$ we need to rule out a pathological situation.

\begin{lemm}\label{lconv}
Under the assumptions of Lemma \ref{lemmaalphaauriga}, let $S:=\Sigma-D^*$ with $D^*$ being the limit of the sequence of diagonal matrices $D_k$ produced by Algorithm 1.
If $S$ has $n$ distinct eigenvalues then
the sequence of rank $r$ matrices $L_k$ produced by Algorithm 1 converges to a rank $r$
matrix $L^*$.
\end{lemm}
\proof
Let $s_1> s_2 > ... > s_n$ be the eigenvalues of $S$ arranged in decreasing order.
By continuity of the eigenvalues, for a sufficiently large $k$, $S_k:=\Sigma-D_k$
has $n$ distinct eigenvalues $s_{k,1}> s_{k,2} > \dots > s_{k,n}$ and
$\lim_{k\rightarrow\infty}s_{k,i}=s_i$.

According to \cite[Chapter 9, Theorem 8]{lax07}, for each $i=1,\dots,n$, we can select an eigenvector
(and hence a normalized eigenvector) $v_{k,i}$ of $S_k$ associated with the eigenvalue
$s_{k,i}$ in such a way that $v_{k,i}$ converges to a normalized eigenvector of $S$ associated with the eigenvalue $s_i$.
Now recall that
$$
L_{k+1}=P_{\mathcal{L}_{n,r}}(S_k)=U_k \diag (f_l (s_{k,1}), \, ..., \, f_l(s_{k,n}) )U_k^\top
$$
where the $i$-th column of $U_k$ is a normalized eigenvector of $S_k$ associated with the eigenvalue $s_{k,i}$.
As a normalized eigenvector is unique up to its sign, we have
$U_k=V_k\Delta_k$ with  $V_k:=[v_{k,1}\mid v_{k,2} \mid\dots\mid v_{k,n}]$
and $\Delta_k$ is a diagonal matrix whose diagonal entries can only be $\pm 1$.
We easily see that  the contribution of the $\Delta_k$  cancels and we have
$$
L_{k+1}=V_k \diag (f_l (s_{k,1}), \, ..., \, f_l(s_{k,n}) )V_k^\top
$$
so that $L_{k+1}$ is given by the product of three matrices each one of which converges as $k$ tends to infinity.
\qed

\begin{propo}\label{propconv}
Assume that the hypothesis of Lemma \ref{lconv} holds and that the matrix $L^*$ defined in the same lemma has rank $r$. Assume also that the tangent space of $\mathcal{L}_{n,r}$ at
$L^*$ does not contain diagonal matrices. Then the sequence ($D_k,L_k$) produced by Algorithm 1 converges to a point corresponding to a  local minimum of the cost function.
\end{propo}
\proof
By the previous results, we know that $D_k$ converges to $D^*$ and $L_k$ converges to $L^*$.
Assume by contradiction that ($D^*,L^*$) is not a  minimum.
Then,  for  any $\varepsilon>0$, there exists $\delta D$ and $\delta L$ such that
$\|\delta D\|_F<\varepsilon$, $\|\delta L\|_F<\varepsilon$, $(L+\delta L)\in {\mathcal{L}_{n,r}}$,
$(D+\delta D)\in {\mathcal{D}_n}$ and
$$
\|\Sigma -L^*-D^* \|_F^2> \|\Sigma -L^*-D^* -\delta L- \delta D \|_F^2.
$$
Now let $\delta T$ be the projection of $\delta L$ on the tangent space of $\mathcal{L}_{n,r}$ at
$L^*$.  For a sufficiently small  $\varepsilon$ we have
$$
\|\Sigma -L^*-D^* \|_F^2\geq \|\Sigma -L^*-D^* -\delta T- \delta D \|_F^2.
$$

By setting $R:=\Sigma -L^*-D^*$ and computing the Frobenius norms in the previous formula, we get
$$
2(\tr[R\delta T] + \tr[R\delta D]) - \|\delta T+ \delta D \|_F^2\geq 0.
$$
By assumption $\delta T+ \delta D\neq 0$ so that
at least one of the two quantities  $2\tr[R\delta L]$ and  $2\tr[R\delta D]$ is positive.
In the case of $\tr[R\delta D]>0$ we have that for all $\kappa$ sufficiently small,
\bea\nn
\min_{D\in\mathcal{D}_n} \Vert {\Sigma} - L^*- D\Vert_F^2 &\leq& \Vert {\Sigma} - L^*- D^*-\kappa \delta D\Vert_F^2\\
\nn
&=& \Vert R\Vert_F^2 + \kappa^2 \Vert \delta D\Vert_F^2- 2\kappa \tr[R\delta D] \\
\nn
&<&\Vert R\Vert_F^2
\eea
which is contradiction because we know that the algorithm converged so that
$\min_{D\in\mathcal{D}_n} \Vert {\Sigma}- L^*- D\Vert_F^2=\Vert R\Vert_F^2$.

In the case of $\tr[R\delta T]>0$ we have that
\beq
\min_{L\in\mathcal{L}_{n,r}} \Vert {\Sigma} - L- D^*\Vert_F^2 \leq  \Vert {\Sigma} - D^* -P_{\mathcal{L}_{n,r}}(L^*+\kappa\delta T )  \Vert_F^2
\eeq
where $P_{\mathcal{L}_{n,r}}(\cdot)$ is the projection onto $\mathcal{L}_{n,r}$.
Thus we have
$P_{\mathcal{L}_{n,r}}(L^*+\kappa\delta T ) =L^*+\kappa\delta T + E$
where $\lim_{\kappa\rightarrow 0} \Vert E\Vert_F/\kappa =0$.

Thus, for $\kappa>0$ sufficiently small, we have
\bea\nn
q&:=& \Vert {\Sigma} - D^* -P_{\mathcal{L}_{n,r}}(L^*+\kappa\delta T )  \Vert_F^2 \\
\nn
&=&\Vert {\Sigma} - D^* -L^* - \kappa\delta T - E  \Vert_F^2
\\
\nn
&=&\Vert R  \Vert_F^2+ \kappa^2 \Vert \delta T  \Vert_F^2 +   \Vert E  \Vert_F^2
 -2 \kappa\tr(R\delta T) -2\tr(R E) +2\kappa\tr(\delta T E)\\
\nn
&<&\Vert R  \Vert_F^2.
\eea
In conclusion, we have
\beq
\min_{L\in\mathcal{L}_{n,r}} \Vert {\Sigma} - L- D^*\Vert_F^2 < \Vert R  \Vert_F^2,
\eeq
that, as in the previous case leads to a contradiction.
\qed

\begin{rem}
We believe that the assumption of Lemma \ref{lconv} can be weakened  that the results still hold
if $s_r > s_{r+1}$ where $s_1\geq s_2 \geq ... \geq s_n$ are the eigenvalues of $S$ repeated according to their algebraic multiplicity and arranged in decreasing order.
The proof of this fact seems, however, very delicate because of some issues on the continuity of eigenspaces under small perturbations.
\end{rem}

\begin{rem}
It is quite intuitive that the conditions of Proposition \ref{propconv} are not very stringent: in fact in all the practical situations that we have studied in simulations those  conditions are satisfied.
\end{rem}

\section{An Alternating Projection Type Algorithm}\label{pojection-approach}

In this section we present our algorithm under a different perspective that may be useful in addressing questions on the properties of the proposed method.
In fact, by suitably translating ${\mathcal{D}_n}$, we  easily see that this method can be viewed as an alternating projection type algorithm for which a very rich literature has been developed.
To this aim,
define
\beq
\begin{aligned}
\mathcal{\tilde{D}}_n:=\Sigma-\mathcal{D}_n=\lbrace X\in\mathbf{Q}_n: \ofd(X)=\ofd({\Sigma}), \, X_{ii}\leq {\Sigma}_{ii}, i=1,...,n \rbrace
\end{aligned}
\eeq
and notice that the projection in this affine set is easily obtained as:
\begin{equation}
P_{\mathcal{\tilde{D}}_n}(X):=
\begin{cases}

X_{ij}=X_{ij} & \text{for } i= j \wedge X_{ii}<{\Sigma}_{ii} \\
X_{ij}={\Sigma}_{ij} & \text{for } (i= j \wedge X_{ii}\geq {\Sigma}_{ii}) \vee i\neq j.
\end{cases}
\end{equation}
We consider now the sequences $L_k$ and $D_k$ produced by our algorithm.
We recall that our  $D_k$  is given by $D_k = P_{\mathcal{D}_n}({\Sigma}-L_k)$.
By taking this formula into account, a direct computation shows that the matrix
$\tilde{D}_k:=P_{\mathcal{\tilde{D}}_n}(L_k)$ may be written as
$\Sigma-D_k$ so that, in view of the formula $L_k= P_{\mathcal{L}_{n,r}}({\Sigma}-D_{k-1})$,
we immediately get that
\[L_{k+1}=P_{\mathcal{L}_{n,r}}(P_{\mathcal{\tilde{D}}_n}(L_k))\]
which shows that the iteration for $L_k$ is the result of  an alternating projection algorithm.
These kind of algorithms burst a long tradition which dates back to Von Neumann in the '30s.
While for alternating projection onto convex sets the convergence results are well established, for the non-convex case much less is known.
In our case $\mathcal{\tilde{D}}_n$ is a convex set of dimension $n$, but the set $\mathcal{L}_{n,r}$ is a non-convex embedded manifold of $\mathbb{R}^{n\times n}$ with dimension $nr-r(r-1)/2$ and it is smooth  at those points for which the rank is exactly $r$.
In \cite{lewis2008alternating} a proof of local convergence (at a linear rate) for alternating projection onto smooth manifolds is provided  under the assumption of transversal intersection.
In our case, transversal intersection cannot hold when $r$ is small with respect to $n$ but it may be possible to generalise that approach to provide a further analysis of the algorithm properties and, in particular, of its convergence rate.

Finally, the set $\tilde{D}_n$ is particularly interesting because of the following interpretation that is particularly  evident when $r$ is such that $\Sigma$ can be decomposed exactly as  $L^*+D^*$ so that $(L^*,D^*)$ is clearly an optimal solution of \eqref{the_problem0}. In this case, $D^*=\Sigma-L^*$ and thus $\Sigma-L^*\in\mathcal D_n$. The latter condition is equivalent to the condition $L^*\in\tilde{\mathcal D_n}$
Therefore the problem \eqref{the_problem0} can reformulated only in terms of $L$ as follows:
\begin{equation}
\begin{aligned}
\label{the_problem}
L^*:=\argmin_{L\in \mathcal{L}_{n,r} \cap \mathcal{\tilde{D}}_n,} \quad & \Vert {\Sigma} -L\Vert_F^2
\end{aligned}
\end{equation}
It is worth noting that the objective function in \eqref{the_problem} does not take into account the covariance matrix of the idiosyncratic noise, i.e. such a matrix is understood as the covariance matrix of a noise random vector. The latter is in the same spirit of  \cite{stock1998diffusion} wherein the factor loading matrix $A$ is given by solving a least squares problem for the linear regression model \eqref{fact_model} and the idiosyncratic component is treated as noise.

\section{Conclusions}
 We have proposed an alternating minimization algorithm for decomposing a covariance matrix as sum of a low rank matrix, whose maximal rank is a priori fixed, plus a diagonal matrix. The latter minimizes the residue among the covariance matrix and the additive decomposition. Simulation results showed that the algorithm performs extremely well and converges very rapidly to the solution. Finally, we have proved that, under reasonable assumptions, such algorithm converges to a solution which is a local minimum for the residue.

\vspace{\fill}\pagebreak


TITLE: Bivariate normal concave downwards
QUESTION [0 upvotes]: The Hessian matrix of the bivariate standard normal distribution (where both standard normal distributions are independent) is positive definite. Yet, it concaves downwards. If the Hessian is positive definite, shouldn't that mean that the function concaves upwards? (Positive definite Hessian implies increasing gradient implies concave upwards.)
Here are the calculations:
Let
\begin{equation}
\phi_2(x, y) = \phi(x)\phi(y) = \frac{1}{2 \pi}e^{-\frac{x^2 + y^2}{2}}
\end{equation}
Observe that $\frac{\partial \phi(x, y)}{\partial x} = -x\phi(x, y)$. The other derivatives follow similar patterns.
Then, the Hessian is
\begin{equation*}
        \begin{matrix}
         \begin{pmatrix}
             \phi_2(x, y)(x^2 - 1) & \phi_2(x, y)xy \\
                \phi_2(x, y)xy & \phi_2(x, y)(y^2 - 1) \\
            \end{pmatrix}
        \end{matrix}
\end{equation*}
which means that the determinant is $\phi_2^2(x, y)(1 - x^2 - y^2)$. Therefore, in the unit circle where $x < 1$ and $y < 1$, the Hessian is positive definite (because the determinant of the Hessian is positive), which is strange because a positive definitive Hessian implies concaving upwards.

REPLY [1 votes]: Your criterion for positive definiteness is incorrect. A positive determinant is necessary but not sufficient for positive definiteness (Sylvester's criterion for positive definiteness of a symmetric matrix requires all leading principal minors to be positive.)
At the point $(0,0)$ the Hessian matrix for the bivariate standard normal is
$
H:=\left(\begin{matrix}-1&0\\0&-1\end{matrix}\right)
$.
This matrix is in fact negative definite. This follows easily from the definition: Given a nonzero $x:=(a,b)^T$, we calculate $
x^THx=-(a^2+b^2)<0$.


TITLE: Proving Sequent Calculus Statement
QUESTION [1 upvotes]: I have to prove the sequent
$$\vdash (\lnot A \lor \lnot B) \to \lnot (A \land B)$$
using the inference rules for natural deduction listed here (pp. 7-8).
I'm super new to natural deduction and sequent calculus. I've attached my work so far below. We need to prove the statement at the bottom, and we're basically working backwards. Honestly, I have no idea

what the next step is, or
how to know when I'm done.

The last step I did was Not Introduction, but that's kinda just because I didn't see any other rule that matched the form of $\Gamma \vdash \lnot P$. Again, complete newbie to this, really don't know what I'm doing or why. I've tried searching things and can understand some simpler formulas (like, just implications and conjunctions/disjunctions), but I don't know how to handle Not.

REPLY [1 votes]: What you did so far goes in the right direction. Now, how to prove the sequent $\lnot A \lor \lnot B, \, A \land B \vdash  \bot $ where you get stuck?
Let us first consider the problem from an informal point of view. Proving the sequent $\lnot A \lor \lnot B, \, A \land B \vdash  \bot $ intuitively means that a contradiction must emerge from the assumptions $\lnot A \lor \lnot B$ and $A \land B$, using the inference rules of natural deduction listed here (pp. 7-8). This is where the rule $\lor_\mathrm{elim}$ plays a crucial role. From the assumption $\lnot A \lor \lnot B$, it follows that there are two cases:

either $\lnot A$ holds, but this is in contradiction with the fact $A$ holds ($A$ is derived from the assumption $A \land B$),

or $\lnot B$ holds, but this is in contradiction with the fact $B$ holds (again, $B$ is derived from the assumption $A \land B$).

In both cases, we have the contradiction we are looking for.
Formally, a rigorous derivation in natural deduction using the inference rules listed here (pp. 7-8) is below.
For the sake of readability, I split the derivation in three parts.
The first one below is the "bottom" part of the derivation. You have to imagine that the derivation $\Pi_A$ is above the sequent $\lnot A \lor \lnot B, \, A \land B, \, \lnot A \vdash \bot$, and that the derivation $\Pi_B$ is above the sequent $\lnot A \lor \lnot B, \, A \land B, \, \lnot B \vdash \bot$.
\begin{align}
  \dfrac
  {\overline{\lnot A \!\lor\! \lnot B, A \!\land\! B \vdash \lnot A \!\lor\! \lnot B}^\mathrm{hyp}
   \
   \overset{\Pi_A}{\lnot A \!\lor\! \lnot B,  A \!\land\! B, \lnot A \vdash \!\bot}
   \ \ \
   \overset{\Pi_B}{\lnot A \!\lor\! \lnot B,  A \!\land\! B, \lnot B \vdash \!\bot}}
  {\dfrac
   {\lnot A \lor \lnot B, \, A \land B \vdash \bot}
   {\dfrac
    {\lnot A \lor \lnot B \vdash \lnot (A \land B)}
    {\vdash (\lnot A \lor \lnot B) \to \lnot (A \land B)}
    \to_\mathrm{intro}
   }
   \lnot_\mathrm{intro}
  }\!
  \lor_\mathrm{elim}
\end{align}
where
\begin{align}
\Pi_A = \
\dfrac
{\dfrac{}{\lnot A \lor \lnot B, \, A \land B, \, \lnot A \vdash \lnot A}^\mathrm{hyp}
 \quad
 \dfrac
 {\overline{\lnot A \lor \lnot B, \, A \land B, \, \lnot A \vdash A \land B}^\mathrm{hyp}}
 {\lnot A \lor \lnot B, \, A \land B, \, \lnot A \vdash A}
 \land_\mathrm{elim}
}
{\lnot A \lor \lnot B, \, A \land B, \, \lnot A \vdash \bot}
\lnot_\mathrm{elim}
\end{align}
and
\begin{align}
\Pi_B = \
\dfrac
{\dfrac{}{\lnot A \lor \lnot B, \, A \land B, \, \lnot B \vdash \lnot B}^\mathrm{hyp}
 \quad
 \dfrac
 {\overline{\lnot A \lor \lnot B, \, A \land B, \, \lnot B \vdash A \land B}^\mathrm{hyp}}
 {\lnot A \lor \lnot B, \, A \land B, \, \lnot B \vdash B}
 \land_\mathrm{elim}
}
{\lnot A \lor \lnot B, \, A \land B, \, \lnot B \vdash \bot}~\lnot_\mathrm{elim}
\end{align}