text:
\begin{document}

\maketitle
\begin{abstract}
In this paper, we describe a low-rank matrix completion method based on matrix decomposition.

An incomplete matrix is decomposed into
submatrices which are filled with a proposed trimming step and then are recombined to form a low-rank completed matrix. The divide-and-conquer approach can significantly reduce computation complexity and storage requirement. Moreover,

the proposed decomposition method can be naturally incorporated into any

existing matrix completion methods to attain further gain.
Unlike most existing approaches, the proposed method is not based on norm minimization nor SVD decomposition. This makes it possible to be applied beyond real domain and can be used in arbitrary fields including finite fields.

\end{abstract}

\section{Introduction}

Consider a large matrix with only a small portion of known entry, an interesting problem is to fill the missing entry assuming the matrix has low-rank. This problem
has several interesting applications including the so-called collaborative filtering problem \cite{srebro2004learning}. An example is the famous Netflix challenge where a huge matrix is used to represent the rating of a movie given by a user. Of course, a typical user will only rate very few movie titles. Therefore, an algorithm will be needed to complete the matrix to predict the ratings of all movies among all users.

It has been shown theoretically that under certain assumptions the matrix can be recovered with very high accuracy \cite{candes2009exact,candes2009power,keshavan2009matrix}. Their approaches convert the rank minimization problem into a nuclear norm minimization problem instead and thus can be solved using semidefinite program (SDP). However, the complexity still grows rather rapidly with the size of the matrix $n$ ($\sim n^3$). Several efficient algorithms have been proposed including Singular Value Thresholding (SVT) \cite{SVT:08}, Atomic Decomposition for Minimum Rank Approximation (ADMiRA) \cite{ADMiRA:09}, Fixed Point Continuation with Approximate (FPCA) \cite{FPCA:09}, Accelerated Proximal Gradient (APG) \cite{APGL:09}, Subspace Evolution and Transfer (SET) \cite{dai2009set}, Singular Value Projection (SVP) \cite{SVP:09}, OptSpace \cite{keshavan2009matrix}, and LMaFit \cite{LMaFit:10}, where OptSpace and SET are based on Grassmann manifold optimization, SVT and SVP uses iterative hard thresholding (IHT) to facilitate matrix shrinkage, FPCA utilizes Bregman iterative algorithm and Monte Carlo approximate SVD, and LMaFit adopts successive over-relaxation (SOR).

In this paper, we propose a decomposition method to allow very efficient divide-and-conquer approach when known entries are relatively very few.
A simple ``trimming'' method is proposed to recover the decomposed ``cluster'' matrix.
However, the decomposition method can also be combine with any other existing matrix completion techniques to yield further gain.
One advantage of the proposed approach is that unlike most existing approaches it does not utilize SVD but only relies on basic vector operations. Therefore, the approach is immediately applicable to matrices of any field (including finite field matrices). This opens up opportunities for new applications.

The rest of the paper is organized as follows. In the next section, we will fix our notation, describe the problem precisely, and present several properties to be used in the later sections. Sections \ref{sect:diag} and \ref{sect:sdiag} will describe the decomposition procedures and present our main results. Section \ref{sect:trimming} will describe the trimming process.

\section{Minimum Rank of Incomplete Matrix}

Let us start with a few notes on our notation.
When things are clear, lines of partition in matrices will not be shown; the $?$ sign may represent an unknown entry, a row or column of unknown entry, a matrix of unknown entry, etc; and similar for the $0$ sign.

Given a finite size matrix $M$ over field $\mathbb{F}$ to be completed, let
\begin{equation}
 S(M)=\{\bar M | \bar M \mbox{ is a completion of } M\}.
\label{rick:1}
\end{equation}

If $ M$ is already completed, then $ S(M)=\{M\}$.
We define
\begin{equation}
mr(M)=\min_{\bar M \in S(M)} \mbox{rank} \bar M.
\label{rick:2}
\end{equation}

Such minimum exists

 because $ \mbox{rank} (S(M)) \subset \mathbb{N}$ and hence $ \exists \bar M \in S(M)$ such that
\begin{equation}
\mbox{rank} \bar M= mr M.
\label{rick:new1.3}
\end{equation}
If $M= \begin{pmatrix} A & B \\ C & D \end{pmatrix}$,
\begin{equation}
 \exists \bar A \in S(A) \mbox{ such that } mr(M) = mr
\begin{pmatrix} \bar A & B \\ C & D \end{pmatrix},
\label{rick:new1.4}
\end{equation}
as we can always find $\bar A$ from $\bar M$ in \eqref{rick:new1.3}. We list other properties about $mr(M)$ that will be quoted:

\begin{align}
 &mr (M) \le \mbox{rank} \bar M,  \qquad \forall \bar M\in S(M)
\label{rick:3} \\
& mr(M) \le mr(P) \mbox{ if  $P$ is any partial completion of  $M$ }
\label{rick:4} \\
 &mr([A | B]) \le mr A + mr B
\label{rick:5}  \\
& mr M^t = mr M
\label{rick:6} \\
 &mr \left( \left[\begin{matrix} A & B \\ C & D \end{matrix}\right]\right) \ge mr(A)
\label{rick:7} \\
&mr M = mr N \mbox{ if $ N$ can be obtained from $ M$ through
interchanging of columns/rows.}
\label{rick:8}
\end{align}

\subsection{Junk Row and Junk Column}
\begin{definition}
A row(column) contains entirely either zero or unknown  will be refered as a junk row(column).
\label{def:9}
\end{definition}

Certainly, we have

\begin{equation}
\label{rick:10}
 mr(J)=0 \mbox { if }  J \mbox{ is a junk row (column), }
 \end{equation}
since we can always complete $J$ entirely by zero entries.

\begin{thm}
\label{thm:1}
Let $ M=[J| N]$ where $ J$ is a junk column, then
$ mr M= mr N$.
\end{thm}

\begin{proof}
By \eqref{rick:5},
we have
$ mr M \le mr J + mr N = mr N$.

On the other hand $ mr M \ge mr N$ by \eqref{rick:7}.
Hence
$ mr M=mr N$.
\end{proof}

Thanks to \eqref{rick:6}, we have the following corollary:

\begin{cor}
\label{cor:1}
$ mr \left( \left[ \begin{matrix} J \\ N \end{matrix} \right] \right) = mr N$ if $ J$ is a junk row.
\end{cor}

\subsection{Equivalence}
We say $ M$ is equivalent to $ N$ and write $ M \sim N$ iff $ N$ can be
obtained from $ M$ through row interchanging, column interchanging,
and junk rows and columns deletion and augmentation.

By Theorem \ref{thm:1} and

and \eqref{rick:8}, we have
\begin{equation}
 M \sim N \Rightarrow mr M = mr N.
\label{rick:11}

\end{equation}

\section{Unknown-diagonalization}
\label{sect:diag}
Define
\begin{equation}
\mbox{u-diag}(B_1,B_2,\cdots,B_n)\triangleq \left[\begin{matrix}
B_1 & \mbox{\bf ?} &\cdots & \mbox{\bf ?}\\
\mbox{\bf ?} &B_2&\cdots&\mbox{\bf ?}\\
\cdots&&\cdots\\
\mbox{\bf ?}&\mbox{\bf ?}&\cdots&B_n\end{matrix}\right]. \label{rick:12}
\end{equation}

We say $M$ is $\mbox{u-diagnonalizable iff } M \sim \mbox{u-diag} (A,B)$,
and both $ A$ and $B$ contain at least one nonzero known entry.

\begin{thm}
Let $ M\sim \mbox{u-diag}(B_1,\cdots,B_n)$, then
$ mr M= max_{1\le i\le n} mr (B_i)$.
\label{thm:2}
\end{thm}

\begin{proof}
By \eqref{rick:11} and induction, all we need to show is when
$ M=\mbox{u-diag}(A,B)$.
Let $ mr A=a$ and $ mr B=b$, by \eqref{rick:7} we have
\begin{equation}
 mr M \ge max (a, b) \label{rick:13}
\end{equation}

Let $ \bar A$, $ \bar B$ be completions of $ A$ and $ B$, respectively
such that
$ \mbox{rank} \bar A= a$ and $ \mbox{rank} \bar B=b$ (c.f. \eqref{rick:new1.3}), then by \eqref{rick:4},
\begin{equation}
mr M\le mr \left( \mbox{u-diag}(\bar A, \bar B) \right)
\label{rick:14}
\end{equation}

Combining \eqref{rick:13} and \eqref{rick:14}, we conclude that
\begin{equation}
 max(a,b) \le mr M \le mr \left( \mbox{u-diag}(\bar A, \bar B) \right).
\label{rick:15}

\end{equation}

By \eqref{rick:8}, we can simply assume the first $ a$ columns
of $ \bar A$ form a basis of $ Col \bar A$ ; and the first $ b$ columns
of $ \bar B$ form a basis of $ Col \bar B$.

Without loss of generality, let us assume $ a\ge b$. We complete the
matrix $ \mbox{u-diag}(\bar A,\bar B)$ by filling up the columns:
\begin{align}
\left[ \begin{matrix}
        \bar A_i \\ \mbox{\bf ?}
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
        \bar A_i\\ \bar B_i
       \end{matrix}
\right],
\nonumber
\\
\left[ \begin{matrix}
        \mbox{\bf ?} \\ \bar B_i
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
        \bar A_i\\ \bar B_i
       \end{matrix}
\right], & \mbox{for $1 \le i \le b$},
\label{rick:new2.5}
 \\
\left[ \begin{matrix}
        \bar A_i \\ \mbox{\bf ?}
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
        \bar A_i\\ \bar 0
       \end{matrix}
\right], & \mbox{for $b+1 \le i \le a$}.
\label{rick:new2.6}
\end{align}
For $i > a$, we make use of the fact that $\bar A_i$ is a linear combination of
$\{ \bar A_k | k \le a \}$. We fill
\begin{align}
\left[ \begin{matrix}
        \bar A_i \\ \mbox{\bf ?}
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
        \bar A_i\\ \sum_{k=1}^a a_{i,k}\bar B_k
       \end{matrix}
\right], & \mbox{for $i > a$}
\label{rick:new2.7}
\end{align}
Similarly,
\begin{align}
\left[ \begin{matrix}
        \mbox{\bf ?} \\ \bar B_i
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
       \sum_{k=1}^b b_{i,k}\bar A_k \\   \bar B_i
       \end{matrix}
\right], & \mbox{for $i > b$},
\label{rick:new2.8}
\end{align}

Now we have a completed $ \mbox{u-diag}(\bar A,\bar B)$ and the first $ a$ of its columns
$ \left[\begin{matrix}\bar A_1\\\bar B_1\end{matrix}\right]$, $ \left[\begin{matrix}\bar A_2\\\bar B_2\end{matrix}\right]$,$ \cdots$,$ \left[\begin{matrix}\bar A_b\\\bar B_b\end{matrix}\right]$,$ \left[\begin{matrix}\bar A_{b+1}\\{\bf 0}\end{matrix}\right]$,$ \cdots$,$ \left[\begin{matrix}\bar A_a\\{\bf 0}\end{matrix}\right]$ form a basis for its column space. Hence it has rank $a$. By \eqref{rick:3}, $ mr \left( \mbox{u-diag}(\bar A,\bar B) \right)\le a$ and hence by \eqref{rick:15} and $ max(a,b)=a$, we get $ mr M=a=mr \left( \mbox{u-diag} (A,B) \right)$ as wanted.
\end{proof}

\begin{rem}
\label{rem:2.1}
Suppose $A$ has been completed by $\bar A$ and $a= \mbox{rank}(\bar A)$. Then if the number of column of $B=n \le a$, then we can complete $B$ arbitrarily and do the completion in \eqref{rick:new2.5}-\eqref{rick:new2.8} as if $b=n$. More generally, given $M=\mbox{u-diag}(\bar A,B)$ with $\bar A$ is completed. Then the completing process of $B$ can be stopped once we know that the final $\mbox{rank} B$ will not be greater than $\mbox{rank} \bar A$ no matter how we do the remaining completion on $B$. For example, if $\mbox{size}(B) \le \mbox{rank} \bar A$, where
\begin{align}
 \mbox{size}(B) = \min (\mbox{number of column of $B$}, \mbox{number of row of $B$}),
\label{rick:new2.9}
\end{align}

then we can complete $B$ arbitrarily to start with.
\end{rem}

\subsection{Percolation and Clusters}

We would call those $B_1,\cdots B_n$ in Theorem \ref{thm:1} as clusters. In other words, clusters are matrices that cannot be u-diagonized. They are not the clusters in the 2-d square lattice, where each point, not counting the edgy one, has 4 neighbors. In our case, each entry in an $n \times m$ matrix has $n + m - 2$ neighbors, from the view of percolation. Despite that difference, the two models share the same percolation threshold at $p \sim 0.6$ \cite{stauffer1994introduction}, where $p$ is the occupation rate. That means if about $60 \%$ of our entries are known, then there is probably one cluster left and the matrix cannot be u-diagonalized. We estimate the number of clusters as the size of the matrix and the number of known entries vary through Monte Carlo simulation and the results are shown in Figure \ref{fig:clusters}. We can see the number of clusters increases as as the number of known entries increases and is peak when the occupation ratio is at about $0.7$ regardness the size of the matrix.

\begin{figure}
\begin{center}
 \includegraphics[width=5in]{num_cluster_vs_K_N10000-crop.pdf}
\caption{The number of clusters versus $k/n$, where $k$ is the number of known entries and the size of the matrix is $n \times n$.}
\label{fig:clusters}
\end{center}
\end{figure}

\subsection{A Decomposition Algorithm}
\label{sect:2.2}

First of all, we set all junk rows and junk columns of the given matrix to zero and blackout them. Now we are working with a junk-free matrix.

We create a row set and a column set for the first cluster. Then put the row position of the first row into the row set; the columns positions of the columns with known entries in that row into the column set. Thanks to \eqref{rick:8}, sorting is not necessary. We black out the row against repeating searches. For each new comers of the column set, we search vertically for its known entries and put the corresponding row position into the row sets. After that we black out the searched columns. Now the row set may have new comers. We enlarge the column set in the way that we enlarged the row set. Both sets keep growing until there is no more new comer.

Then we create another column set and another row set. Repeat the procedure for the next cluster, if the remaining matrix has not been blackout to void.

\section{Sub Unknown-diagonalization}
\label{sect:sdiag}
\begin{definition}
\label{definition:3.1}
Let $u$ be a vector and $A$ be a matrix. We define
 $zero(u, A)$ as $1$ if  $\forall [\bar u | \bar A]\in S([u | A])$
 (c.f. \eqref{rick:1})
with $\mbox{rank} [\bar u | \bar A]= mr [u | A]$, then $\bar u= {\bf 0}$.
Otherwise, we define $zero(u, A)$ as $0$.
\end{definition}

\noindent
E.g.
$
zero \left(
\begin{pmatrix} 0 \\ {\bf ?} \end{pmatrix},
\begin{pmatrix} 1 \\ 1 \end{pmatrix}
\right)=1$,
$zero\left(
\begin{pmatrix}0 \\ {\bf ?} \end{pmatrix},
\begin{pmatrix}0 \\ 1 \end{pmatrix}
\right)=0$.
{\color{black} Notice that the necessary condition (not sufficient) for $zero (u, A)=1$ is
that $u$ must be a junk column.}
Hence we have
$zero(u,A) \le junk(u)$, where $junk(u)=\left\{\begin{array}{l l}
1, &               \mbox{if $u$ is junk,} \\
0, & otherwise.
                                        \end{array}\right.$.

\begin{thm}
\label{thm:5}
Let $M=
\begin{pmatrix}
v_1 & B1& &        \mbox{\bf ?} \\
v_2  & &     B2  \\
\cdots & & &           \cdots \\
v_n &     \mbox{\bf ?} &&          Bn
\end{pmatrix}
$, where $B_i$ are matrices, $v_i$ are vectors that
$\begin{pmatrix}
v_1 \\
v_2 \\
\cdots \\
v_n
 \end{pmatrix}
$ is not a junk column, then
$mr M= \max_i (mr[v_i | B_i] + zero (v_i, Bi))$.
\end{thm}

\begin{proof}
We will show the case when $n=2$; cases of higher $n$ are easy to generalized. Let
$M=\begin{pmatrix} u & A & \mbox{\bf ?} \\
v &\mbox{\bf ?} & B\end{pmatrix}$ with
$\begin{pmatrix}
  u \\v
 \end{pmatrix}
$ is not a junk column.
If $zero (u, A)=0$, then we have
\begin{align}
mr M & \ge mr [u | A]  \nonumber \\
& =mr [u | A] + zero (u, A), \label{rick:25}
\end{align}
where the first inequality is by \eqref{rick:7}.
If $zero (u, A)=1$, then $u$ must be a junk column. Since
$\begin{pmatrix}
u \\
v
 \end{pmatrix}
$ is not a junk column, $v$ must contains a nonzero known
entry $\lambda$.
Let $M' = \begin{pmatrix}
u   &         A\\
\lambda  & {\bf ?} \cdots {\bf ?}
          \end{pmatrix}
$. By \eqref{rick:5} and \eqref{rick:7},
\begin{equation}
 1+mr A \ge mr M' \ge mr A. \label{rick:26}
\end{equation}
So there are only two possibilies for $mr M'$.
Assume $mr M'= mr A$. By
\eqref{rick:new1.3}

we can
pick $[\bar u | \bar A]\in S([u | A])$ such that
\begin{equation}
mr A= mr M' = mr \begin{pmatrix} \bar u   &   \bar A \\
\lambda   & {\bf ?}\cdots {\bf ?}\end{pmatrix}
\ge
mr [ \bar u | \bar A ].
\label{rick:27}
\end{equation}
Therefore,
$mr A =mr [\bar u \bar A]$ as
$rank[\bar u \bar A]\ge rank \bar A \ge mr A$. Since $zero (u,A)=1$,

we have
$\bar u={\bf 0}$ (c.f. def \ref{definition:3.1}).
But that will make
$mr M’ = rank \bar A +1$ because $\lambda \neq 0$.
We conclude that $mr M'$ cannot equal to $mr A$ and hence we must have
\begin{align}
mr M & \ge mr M’= mr A+1  &\mbox{(c.f. \eqref{rick:26})} \nonumber \\
     & =  mr [u | A] + zero (u, A). & \mbox{(c.f. Theorem \ref{thm:1})}
\label{rick:new3.4}
\end{align}

Combine \eqref{rick:25} and \eqref{rick:new3.4} and the symmetry between $(u, A)$
and $(v, B)$ we get
\begin{equation}
mr M \ge \max (mr[u | A] +zero (u, A), mr[v | B] + zero (v, B)).
\label{rick:new3.5}
\end{equation}

Without loss of generality, we may assume
\begin{align}
a=mr[u| A] + zero (u,A) \ge mr [v| B] + zero (v, B)=b.
\label{rick:new3.6}
\end{align}

Pick $[\bar A_0| \bar A]\in S([u| A]), [\bar B_0| \bar B]\in S([v| B])$
s.t.
\begin{align}
rank  [\bar A_0| \bar A]& = mr [u | A], \nonumber \\
 rank  [\bar B_0| \bar B]& = mr [v| B].  \qquad \mbox{(c.f. \eqref{rick:2})}
\label{rick:new3.7}
\end{align}

and
\begin{align}
\bar A_0 &\neq  \vec 0,  & \mbox{ if $zero (u, A)=0$},& \nonumber \\
  \bar B_0 &\neq \vec 0, & \mbox{ if $zero (v, B)=0$}.
&\qquad  \mbox{(c.f. def \ref{definition:3.1})}

\label{rick:new3.8}
 \end{align}

Thanks to \eqref{rick:8}, we can assume
$$Span \{\bar A_0, \bar A_1,....,\bar A_{a-1}\}=Col [\bar A_0| \bar A]$$,
$$Span \{\bar B_0, \bar B_1,....,\bar B_{b-1}\}=Col [\bar B_0| \bar B]$$.

Then we complete the rest by filling up the columns:
\begin{align}
\left[ \begin{matrix}
        \bar A_i \\ \mbox{\bf ?}
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
        \bar A_i\\ \bar B_i
       \end{matrix}
\right],
\nonumber
\\
\left[ \begin{matrix}
        \mbox{\bf ?} \\ \bar B_i
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
        \bar A_i\\ \bar B_i
       \end{matrix}
\right], & \mbox{for $1 \le i < b$},
\label{rick:new3.9}
 \\
\left[ \begin{matrix}
        \bar A_i \\ \mbox{\bf ?}
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
        \bar A_i\\ \bar 0
       \end{matrix}
\right], & \mbox{for $b \le i < a$}.
\label{rick:new3.10}
\end{align}
For $i \ge a$, we can pick $\{a_{i,k} \} \in \mathbb{F}$ s.t.
$\bar A_i = \sum_{0 \le k < a} a_{i,k} \bar A_k$ and fill the columns
\begin{align}
\left[ \begin{matrix}
        \bar A_i \\ \mbox{\bf ?}
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
        \bar A_i\\ \sum_{0 \le k < a} a_{i,k}\bar B_k
       \end{matrix}
\right], &
\label{rick:new3.11}
\end{align}
Similarly, for $i \ge b$, we pick $\{ b_{i,k} \} \in \mathbb{F}$ s.t.
$\bar B_i = \sum_{0 \le k < b} b_{i,k} \bar B_k$ and fill the columns
\begin{align}
\left[ \begin{matrix}
        \mbox{\bf ?} \\ \bar B_i
       \end{matrix}
\right] & \Rightarrow
\left[ \begin{matrix}
       \sum_{0 \le k < b} b_{i,k}\bar A_k \\   \bar B_i
       \end{matrix}
\right], &
\label{rick:new3.12}
\end{align}
Then the completed matrix is rank $a$ with the first $a$ columns form a basis of its column space. By \eqref{rick:3},
\begin{align*}
      mr M & \le  a  &\\
          &= rank [\bar A_0 \bar A] + zero (u, A)  &
\mbox{(c.f. \eqref{rick:new3.6}, \eqref{rick:new3.7})} \\
           &=\max (mr[u A] +zero (u, A), mr[v B] + zero (v, B)) &
\end{align*}
as assumed. Together with \eqref{rick:new3.5}, we get
$mr M=\max(mr[u A] +zero (u, A), mr[v B] + zero (v, B))$.

\end{proof}

\begin{rem}
\label{rem:3.1}
 Suppose $[u| A]$ has been completed by $[\bar u| \bar A]$ and we have $a=zero (\bar u, \bar A) + rank[\bar u \bar A]$. Then if the number of column of $[v | B]=n \le a$, we can complete $[v| B]$ arbitrarily and do the completion in \eqref{rick:new3.9}-\eqref{rick:new3.12}
 as if $b=n$. More generally, the completing process of $B$ can be stopped once we know that the final $zero(\bar v, \bar B) + rank [\bar v \bar B]$ won't be greater than $a$ no matter how we do the remaining completion on $[v B]$. For example if $junk(v)+size([v B]) \le a$ (c.f. \eqref{rick:new2.9}, def \ref{definition:3.1}), then we can complete $[v| B]$ arbitarily at the beginning and do the completion \eqref{rick:new3.9}-\eqref{rick:new3.12}.
\end{rem}

\subsection{How to decompose sub u-diagonalizable matrix}

\begin{definition}
 A matrix $C$, not u-diagonalizable, becomes u-diagonalizable after deleting a row or a column is called sub u-diagonalizable. The row (column) is called conjoined row (column).
\end{definition}

For example,
$\begin{pmatrix}
  v_1 \\ v_2 \\ \cdots \\v_n
 \end{pmatrix}
$
in Theorem \ref{thm:5} is a conjoined column and $M$ becomes
u-diagonalizable without it.

 \begin{definition}
\label{definition:3.3}
Given two vectors $v$ and $w$ of same length, we say $v$ is a {\em donor} for $w$ $(v \succeq w)$ iff all of the unknown positions of $v$ are also unknown in $w$. In order words, after some row interchanging $[w | v]
=
\begin{pmatrix}
 \bar r & \bar d \\
\vec ?  &  n
\end{pmatrix}
$ with $\bar r$ and $\bar d$ are completed. Clearly, $v \succeq  w$ and $w \succeq u$ imply $v \succeq u$. However, if $v \succeq w$ and $w \succeq v$, we do not have $v=w$. Vectors $v$ and $u$ are said to be {\em comparable} if either $v \succeq w$ or $w \succeq v$.
 \end{definition}

\begin{thm}
\label{thm:3.2}
 Conjoined row (column) does not have donors among other rows of the sub u-diag matrix.
\end{thm}

\begin{proof}
 Let $C$ be the sub u-diag  matrix, then it must have the following structure (after some row and column interchanging):
$
C = \begin{pmatrix}
 A & ? \\
           ? & B  \\
           u & v
\end{pmatrix},
$
where $u$ cannot be entirely unknown, otherwise $C$ is u-diagonalizable. Now, rows in $[? | B]$ cannot be donors of $[u | v]$, the conjoined row. Similarly $v$ cannot be entirely unknown and hence, rows in $[A | ?]$ cannot be donors of $[u| v]$ neither.
\end{proof}

Therefore if $C$ is a sub unknown-diagonalizable, we will not miss the chance of decomposing it if we have tested every row and column that does not have a donor. That is to blackout the suspicious row (column) and then carrying out the decomposition mentioned in Section \ref{sect:2.2}. We would like to call the decomposed components as {\em sub-clusters}. For example,
$\begin{pmatrix}
 A \\u
\end{pmatrix}$
and $\begin{pmatrix}
     B \\v
    \end{pmatrix}$
are sub-clusters of the $C$ in the proposition.

Unlike cluster that cannot be further unknown-diagonalized, sub-clusters can be sub unknown-diagonalizable. For example
$C=\begin{pmatrix}
A & ? & ?  \\
 u & v & ?   \\
? & B & ?    \\
? & x & y     \\
? & ? & D
 \end{pmatrix},$
where both $[u | v | ?]$ and $[?| x| y]$ are conjoined rows. In that case, we may first decompose $C$ into sub-clusters
$\begin{pmatrix}
  A \\u
 \end{pmatrix}$
and
$\begin{pmatrix}
  v & ? \\ B & ? \\x & y \\ ? & D
 \end{pmatrix}$.
 Then we may further decompose the later into
$\begin{pmatrix}
  v \\ B \\ x
 \end{pmatrix}$
and
$\begin{pmatrix}
  y \\ D
 \end{pmatrix}$, if necessary.

\section{Trimming}

\label{sect:trimming}
\begin{lem}
 $mr [v | M]=mr M$ if $\forall \bar M \in S(M)$,
$Col (\bar M) \cap S(v) \neq \mathbb{\phi}$,.
\label{lem:4.1}
\end{lem}
\begin{proof}
 By (1.9) we already have
$mr[v| M] \ge mr M$.
Let $\bar M \in S(M)$ such that $rank \bar M=mr M$ (c.f. \eqref{rick:2}).
Then pick a $\bar v \in S(v) \cap Col(\bar M)$. From \eqref{rick:3}, we get
$mr [v | M] \le rank [\bar v | \bar M]=rank \bar M=mr M$.
\end{proof}

\begin{thm}
\label{thm:4}
Let $M_{d_i}$ be the $d_i$-th column of $M$ and donor (c.f. def \ref{definition:3.3})
of a vector $v$ for $1\le i \le t$, such that after some row interchanging,
$[v | M_{d_1} M_{d_2} \cdots M_{d_t}]=
\begin{pmatrix}
\bar r | \bar D \\
                             ?  | N
\end{pmatrix}$
 with $\bar r$ and $\bar D$ are completed.
If $\bar r\in Col (\bar D)$, then $mr[v | M]= mr M$.
\end{thm}

\begin{proof}
Thanks to \eqref{rick:8}, we may start with
$
[v |M]=
\begin{pmatrix}
 \bar r |  \bar D A \\
? | N B
\end{pmatrix}
$

Pick $a_i\in \mathbb F$ such that
\begin{equation}
\bar r=\sum a_i \bar D_i.
\label{rick:new4.1}
\end{equation}
Then $\forall
\bar M=
\begin{pmatrix}
 \bar D \bar A \\
\bar N \bar B
\end{pmatrix}\in S(M)$,

we complete $v$ to
\begin{align}
 \bar v=
\begin{pmatrix}
\bar r \\
        \sum a_i\bar N_i
\end{pmatrix}
\in S(v)\cap Col (\bar M).
\label{rick:new4.2}
\end{align}
(Note that $rank[\bar v | \bar M]= rank \bar M$.)
Now Lemma \ref{lem:4.1} implies $mr[v | M]=mr M$.
 \end{proof}

\subsection{Trimming Process}

We test column by column to see if we can make use of Theorem \ref{thm:4} to trim away some columns from a given matrix, which is probably a sub-cluster mentioned in the previous section. We call this process as {\em column trimming}.
When we find a column satisfying the condition of Theorem \ref{thm:4}, we will mark down the dependency relation between it and its donor (i.e.
\eqref{rick:new4.1}) in order. Then we black it out and go for the next column.

Similarly, we have row trimming. An uninterrupted (c.f. Remarks
\ref{rem:2.1} and \ref{rem:3.1}) trimming process starts with a column trimming followed by a row trimming, or the other way round. Then we carry out these two kinds of trimming one after the other, until there is no more reduction in the matrix.
After the trimmed matrix gets completed, we restore, in reverse order, the blackouts with the completed forms given by \eqref{rick:new4.2}.

The following proposition is interesting in its own right and may be useful for our future study.

\begin{lem}
\label{lem:4.2}
 If $v$ is a column that $S(v)\cap Col(\bar M)=\mathbb{\phi}$, $\forall
\bar M\in S(M)$. Then $mr[v  | M]=mr M+1$.
\end{lem}

\begin{proof}
By \eqref{rick:5} $mr [v| M] \le mr M +1$. Let us complete $[v | M]$ by $[\bar v | \bar M]$ such that
$rank [\bar v | \bar M]= mr [v | M]$ (c.f. \eqref{rick:2}). Since $\bar v \notin Col (\bar M)$, we must have $rank [\bar v | \bar M]= rank \bar M+1$, which implies
$mr [v | M] = rank \bar M + 1 \ge mr M + 1$ (by \eqref{rick:3}).
\end{proof}

\noindent
{\bf Proposition.} {\it
Suppose $M$ is a matrix that every two columns $M_i$ and $M_j$ of $M$ are comparable (c.f. def \ref{definition:3.3}). Then one round of column trimming followed by arbitrarily completion and proper restoration (i.e. \eqref{rick:new4.2}) of the trimmed columns complete $M$ to its minimum rank.  }

\begin{proof}
Let $T$ be the trimmed $M$. Every two columns of $T$ are also columns of $M$ and hence comparable. So there exists a column $T_j$ in $T$ such that $T_j \preceq T_i$ for all $i$. After some row interchanging and column interchanging, we have
$T_j=\begin{pmatrix}
  \bar r \\?
 \end{pmatrix}$
and
$T = \begin{pmatrix}
      \bar r & \bar D  \\
? & N
     \end{pmatrix}$
with $\bar r$ and $\bar D$ are completed. Notice that $\bar r\notin Col (\bar D)$, otherwise $T_j$ has been trimmed away already. Now Lemma \ref{lem:4.2} and \eqref{rick:5} give $mr N=1+ mr
\begin{pmatrix}
 \bar D \\ N
\end{pmatrix}
$. Repeating the argument, we get
$mr N= \mbox{number of columns of }N \ge rank (\bar N)$, $\forall \bar N\in S(N)$. Together with Theorem
\ref{thm:4} and
\eqref{rick:3}, we conclude that $mr M=mr N= rank (\bar N)$, $ \forall N\in S(N)$. Finally, the proper restoration (c.f. \eqref{rick:new4.2}) restores $\bar N$ to $\bar M\in S(M)$ such that $rank (\bar M)=rank (\bar N)=mr (M)$.
 \end{proof}

\subsection{Trimming Process with Approximation}
The trimming process stops when there are no more columns or rows fulfill the condition of Theorem \ref{thm:4}. But we can always make an approximation by blacking out a column or a row as if it fulfill the condition and continue the trimming. With this approximation, the process stops when there is no more unknown left in the trimming matrix. (We may even choose the row or column that has no donor (c.f. Thm \ref{thm:3.2}) to black out and check for u-diagonalization. )

Then we restore the blackouts in reverse order. When we meet a blackout without dependency relation (i.e. \eqref{rick:new4.1}) to restore, we check for the condition of Theorem \ref{thm:4}, again. The first time was with the uncompleted trimming matrix; this time is with the completed restoring matrix. If the condition is fulfilled, we restore the blackout with completed form given by
\eqref{rick:new4.2}. This will not compromise (further) the minimum rank that we can reach. Otherwise we restore the blackout with arbitrary completed form, which may cause one (more) rank deviation from the possible minimum.

\bibliographystyle{IEEEtran}

\bibliography{../../ref}
\end{document}


TITLE: Irreducible non-Abelian subgroup of $\mathrm{U}_n(\mathbb{C})$, containing diagonal matrices
QUESTION [1 upvotes]: Consider an irreducible non-Abelian subgroup $\mathrm{H}$ of group of unitary matrices $\mathrm{U}_n(\mathbb{C})$, that contains the subgroup of diagonal matrices. Does there exist any result regarding the properties or structures of $\mathrm{H}$? The final motivations is to answer this question.
P.S. The question has been updated after some useful comments.

REPLY [3 votes]: 1) First, let $H$ be a closed connected subgroup with this property. Let $D$ the diagonal group in $U(n)$; denote Lie algebras with Gothic letters. Then $$\mathfrak{u}(n)=\mathfrak{d}\oplus \bigoplus_{j<k}\mathbf{R}i(E_{jk}+E_{kj})\oplus\mathbf{R}(E_{jk}-E_{kj}).$$
Let $e_j:D\to\mathbf{C}^*$, $d\mapsto d_j$ be the projection (valued in the unit circle). As $D$-module, the above decomposition of $\mathfrak{u}(n)$ is invariant, $\mathfrak{d}$ has weight $e_j$ and $\mathbf{R}i(E_{jk}+E_{kj})\oplus\mathbf{R}(E_{jk}-E_{kj})$ is $\mathbf{R}$-irreducible with 2-dimensional complexification of weights $\pm e_j-e_k$. Since these are distinct when the pair $(j,k)$ with $j<k$ varies, it follows that any $D$-submodule $M$ of $\mathfrak{u}(n)$ containing $\mathfrak{d}$ has the form
$$M=\mathfrak{d}\oplus \bigoplus_{j<k;(j,k)\in W}\mathbf{R}i(E_{jk}+E_{kj})\oplus\mathbf{R}(E_{jk}-E_{kj})$$
for some subset $W$ of the set of pairs $(j,k)$ with $j<k$. Let $W'$ be the set of pairs $(j,k)$ such that $(j,k)$ or $(k,j)$ belongs to $W$. So $W'$ is symmetric and $$M=\mathfrak{d}\oplus \sum_{(j,k)\in W'}\mathbf{R}i(E_{jk}+E_{kj})\oplus\mathbf{R}(E_{jk}-E_{kj})$$
The condition that $M$ is a Lie subalgebra easily implies that $(j,k),(k,\ell)\in W'$ imply $(j,\ell)\in W'$. Hence $W''$, the union of $W'$ and the diagonal, is an equivalence relation on $\{1,\dot,n\}$. Conversely, for every equivalence relation $W''$ on $\{1,\dots,n\}$,
$$\mathfrak{h}_{W''}=\mathfrak{d}\oplus \sum_{j\neq k;(j,k)\in W''}\mathbf{R}i(E_{jk}+E_{kj})\oplus\mathbf{R}(E_{jk}-E_{kj})$$
is a Lie subalgebra containing $\mathfrak{d}$. The corresponding group is thus the group of block-diagonal matrices with respect to some partition (possibly permuting indices to make it block-wise).
2) Now let $H$ be a closed subgroup containing $D$, possibly not connected. Then $H^0$ has the preceding form, and $H$ normalizes $H^0$. One can check that the normalizer of $H^0$ is always of finite index over $H^0$: indeed, the blocks are precisely the irreducible components of the $H^0$-actions, and are pairwise non-isomorphic $H^0$-modules, hence they are permuted by $H$. That is, this normalizer is the stabilizer of some direct sum according to some partition of indices, possibly permuting blocks.
3) If one wants irreducibility (as I said in a comment, it just complicates the discussion while restricting the scope): it corresponds to the case where $H/H^0$ acts transitively on the set of blocks (this is possible only if all blocks have the same size)
4) The remaining step is to show that any subgroup $H$ containing $D$ is automatically closed. To start with, the connected component of the closure of $H$ is the component-wise stabilizer of some partition $P$ of $\{1,\dots,n\}$.
Let $x=(x_1,\dots,x_m)$ be an $m$-tuple of $H$. Consider the map $D^m\to U(n)$ mapping $(d_1,\dots,d_m)$ to $\prod x_id_ix_i^{-1}$. Let $r_x$ be its rank (maximum rank of its differential over $D^m$). So for some $y=(y_1,\dots,y_m)$, its rank at $y$ is $r_x$. Hence for $x'=(x_1,\dots,x_m,x_m,\dots,x_1)$, its rank at $(y_1,\dots,y_m,y_m^{-1},\dots,y_1^{-1})$ is $\ge r_x$ and moreover the value is $1$. Thus, we can assume that $x$ is chosen such that $r_x$ is maximal and achieved at a point $(y_1,\dots y_m)$ with value $1$. From the maximality it follows that the tangent image is a Lie subalgebra $\mathfrak{l}$ of $\mathfrak{u}(n)$, and it does not depend on the choice of $x$, and the corresponding immersed Lie subgroup $L$ is contained in $H$ and contains $D$. The preceding results concerns Lie subalgebras containing $\mathfrak{d}$ applies, so $\mathfrak{l}$ is
the stabilizer of some partition $Q$ of $\{1,\dots,n\}$ (with $Q\subset P$ since $L\subset \bar{H}^0$). But it it easy to see that if $P\neq Q$ then $\mathfrak{h}_Q$ is not normalized by $\mathfrak{h}_P$. So $P=Q$. Hence $H\supset L=\bar{H}^0$. It follows that $H$ is closed.


TITLE: What exactly is the equivalence class of the product of ideals?
QUESTION [0 upvotes]: I was reading So What Is Class Number 2? by Scott T. Chapman and it said that for two ideals of $R$, $I$ and $J$, we can say that $I\sim J$ if there is $\alpha,\beta\in R$ such that $(\alpha)I=(\beta)J$. Furthermore, we can say $[I]$ is the equivalence class of $I$, $[J]$ is the equivalence class of $J$, and $IJ=\{\Sigma_{i=1}^na_ib_j\ |\ a_i\in I\ \text{and}\ b_j\in J\}$.
Here's where I'm getting tripped up: the paper goes on to claim that $[I]+[J]=[IJ]$.
What does $[IJ]$ actually mean, though?
I think I have a good idea what $[I]$ implies. If, for example, $I$ is the set of even numbers in $\mathbb{Z}$, then $[I]$ would be the set of all ideals $J$ such that $(\alpha)2\mathbb{Z}=(\beta)J$ for integer $\alpha$ and $\beta$. But, to the best of my knowledge, $IJ$ does not appear to be a set. It is a sum. If $I=\langle x,y\rangle$ and $J=\langle z \rangle$, then wouldn't $IJ=xz+yz$? If so, then how can you define an equivalence relation and thus an equivalence class with this? Or am I misunderstanding the definition? Help would be appreciated. Thank you!

REPLY [1 votes]: It appears your confusion is really about what the product of two ideals actually is. You seem to have misunderstood the "set builder" notation as denoting a single element (as per the comments). I note that you also have an undefined $j$ in the definition of $IJ$ in the first paragraph; they should both be $i$.
The notation
$$\left\{ \sum_{i=1}^n a_ib_i\,\Bigm|\, a_i\in I,b_i\in J\right\}$$
does not represent a single element of the given form, but rather all elements that can be written this way. This is (somewhat abusive) shorthand for
$$\left\{ x\,\Bigm|\,\text{there exist a positive integer }n\text{ and elements }a_i\in I,b_j\in J\text{ such that }x=\sum_{i=1}^n a_ib_i\right\}.$$
Thus, if $I=\langle x\rangle$ and $J=\langle y,z\rangle$ are ideals in (the commutative ring with unity) $R$, then $I$ consists of all elements of the form $ax$ with $a\in R$, and $J$ consists of all elements of the form $by+cz$ with $b,c\in R$. The product $IJ$ will consist of all sums of elements of the form $a_ix(b_iy+c_iz)$.
In general, the set of products $ab$ with $a\in I$ and $b\in B$ need not be an ideal, because it may not be closed under sums. For example, if $R=\mathbb{Z}[x,y,z,w]$ with $I=(x,y)$ and $J=(z,w)$, then $xz$ and $yw$ are both products of elements of $I$ with elements of $J$; but $xz+yw$ cannot be expressed as a product $(ax+by)(cz+dw)$ with $a,b,c,d\in\mathbb{Z}$ of an element of $I$ with an element of $J$; hence the need for taking sums of products when defining the product of the ideals.


TITLE: Rudin exercise 3.8 solution
QUESTION [0 upvotes]: Problem . if the sequence $\sum\limits_{}{a_{n}}$ converges and $b_{n}$ is a monotonic bounded sequence then $\sum\limits_{} {a_{n}b_{n}}$ converges.
Solution. Let's use the identity:
$\sum\limits_{n=p}^q {a_{n}b_{n}}$ = $\sum\limits_{n=p}^{q-1} {s_{n}(b_{n}-b_{n-1})}+s_{q}b_{q}-s_{p-1}b_{p}$, where $s_{n}=\sum\limits_{k=1}^n {a_{n}}$. Then $|\sum\limits_{n=p}^{q-1} {s_{n}(b_{n}-b_{n-1})}|<=|\sum\limits_{n=p}^{q-1} {s_{n}}|B$, where B is the difference of upper and lower bound of the sequence $b_{n}$. since $a_{n}$ is convergent we can make $|\sum\limits_{n=p}^{q-1} {s_{n}}|$ as small as we want and therefore we can make this $|\sum\limits_{n=p}^{q-1} {s_{n}(b_{n}-b_{n-1})}|$ as small as we want like $x/3$. Then let's take the absolute value of another term of the right side of the identity.
$|s_{q}b_{q}|<=|s_{q}||b_{up}|$, where $b_{up}$ is the upper bound of the sequence and the same argument here too. we can make this term as small as we want. The same goes to the absolute value of the third term on the right side of the identity. Finally $|\sum\limits_{n=p}^{q} {a_{n}b_{n}}|<=|\sum\limits_{n=p}^{q-1} {s_{n}(b_{n}-b_{n-1})}|+|s_{q}b_{q}|+|s_{p-1}b_{p}|$. We can find $p$ and $q$ to make the right side small enough, completing the proof. Is this right? I am not satisfied what I have written here, I should have stated better in my opinion but the general idea of the proof is shown I suppose.

REPLY [3 votes]: Assume that $\{b_n\}$ is increasing. If $s_n=\sum_{k=1}^n a_k$, then
$$
\sum_{k=1}^n a_kb_k=\sum_{k=1}^n (s_k-s_{k-1})b_k=\sum_{k=1}^n s_kb_k-\sum_{k=2}^n s_{k-1}b_k=s_nb_n+
\sum_{k=2}^{n}s_{k-1}(b_{k-1}-b_{k}).
$$
Now $\{s_nb_n\}$ is convergent, as a product of convergent sequences, while the series
$\sum s_{n-1}(b_{n-1}-b_{n})$
is dominated by the series $\sum M(b_{n}-b_{n-1})$, where $M=\sup_{n\in\mathbb N}|s_n|$, and since $\{b_{n}-b_{n-1}\}$ is absolutely summable,
so is  $\sum s_{n-1}(b_{n-1}-b_{n})$.


TITLE: What sort of operations can be applied on a Hilbert spaces?
QUESTION [1 upvotes]: I was reading the paper No Universal Flipper for Quantum States. In this paper they have tried to prove by contradiction that a universal flipping machine cannot exist. By flipping I mean if I have a qubit $\alpha|0 \rangle +\beta|1\rangle$, then the flipped qubit will be $\beta^*|0 \rangle -\alpha^*|1\rangle$. In the paper they take an entangled system between two parties Alice and Bob, and apply an operation on Bob's side. But they say
that the operation is not necessarily unitary ( its just linear ) but it does preserve the trace of Alice's reduced density operator before and after the operation.
But I thought on closed quantum mechanical systems only unitary operations could be applied. Am I missing something?

REPLY [3 votes]: Correct, the evolution of a quantum system does not necessarily need to be unitary. However, if the system is closed, it needs to be unitary. The latter is due to the fact that we assume the system evolves due to the Schrödinger equation with a self-adjoint Hamiltonian, which induces a unitary evolution operator via Stone's theorem.
For non-closed systems, the idea is the following: Let's suppose you have your system (denote it by a state $\rho$) and maybe some other system, the environment (denoted by a state $\rho_E$). Then you start out with the state $\rho\otimes \rho_E$ and you can perform any unitary operation you want, i.e. $\tilde{\rho}:=U(\rho\otimes \rho_E)U^*$. Now so far, we have dealt with what you learn in your first quantum theory course. The idea is that in our case, we are simply not interested in (or can't keep track of) the evolution of the environment. Hence the state we have at the end is the reduced density matrix of $\tilde{\rho}$ restricted to our system, i.e:
$$ \rho_{final}=\operatorname{tr}_E(U(\rho\otimes \rho_E)U^*) $$
with the partial trace over the environment. Now it turns out that this is (more or less) the most general form of what is called a completely positive map and hence we usually say that those are the maps governing quantum mechanics in closed as well as open systems. In order to keep the normalization of a state, these need to be trace-preserving, hence we say that a general quantum evolution is governed by a completely positive trace-preserving map (CPTP) or quantum channel.
Note that the only point to argue is about $\rho\otimes \rho_E$: It is not so clear that we should have to start with an uncorrelated systems, but there are good reasons to do so.
Now, maybe you don't even know how the environment looks like. These systems are called "open quantum systems" and one way to describe their evolution is by as master equation as in statistical equation (for quantum mechanics, the corresponding type is also often known as Lindblad-equation). However, the corresponding map on density matrices is also a quantum channel and thus of the form outlined above.
Another characterization of quantum channels might be helpful for your understanding: They can be axiomatized (in finite dimensions for simplicity) via:

A quantum channel is a linear map from matrices to matrices (we need linearity physically; see below).
It is a positive map, i.e. it sends positive semidefinite matrices to positive semidefinite matrices (states are positive semidefinite and our map better send states to states).
Moreover, if $T$ is the map, then also $T\otimes \operatorname{id}_n$, where $\operatorname{id}_n$ is the identity is also a positive map for all $n$ (this is what is meant by "complete" positivity: If our map acts just on a subsystem, the natural extension to the whole system should definitely be a valid map).
It preserves the trace of the matrices.

Except for maybe the third property, it should be immediately clear that the other properties should be true for any quantum evolution starting of any system - and these should also be the properties the authors use for their proof.
Linearity: Finally, let's talk about antiunitary evolutions: It turns out that a map as proposed is not linear (it's antilinear), so we have to exclude antiunitary maps, too, because they would violate the linearity of quantum mechanics. This linearity however is established by many experiments and lies at the heart of quantum theory - there is no good reason why we should abandon it.


TITLE: What is the Cayley graph of $(\Bbb Z/2\Bbb Z)\times(\Bbb Z/2\Bbb Z)$?
QUESTION [2 upvotes]: I get that the presentation of the new group, with respect to two generators, would be $(x,y \;|\; x^2= y^2=1)$ but I'm not sure how the actual graph would look. Would it consist of an infinite sequence of loops, where each loop has one edge for the $x$ generator and one for the $y$? But then how would the $x^2$ relations come into play? Any help in my understanding would be greatly appreciated. Thanks!

REPLY [2 votes]: No it's not infinite. You have the four points $1, x, y, xy$ and with $S=\{x,y\}$ as your generating set, you get edges $\{1,x\}, \{1,y\}, \{x,xy\}, \{y, xy\}$.
Also your presentation is incorrect, what you have written is $\Bbb Z/2\Bbb Z * \Bbb Z/2\Bbb Z$, you need to add the relation $xyx^{-1}y^{-1}=1$, otherwise there are infinitely many nodes.
Note each node has degree $2$, which is exactly right since $|S|=2$. So the graph looks like a square with $1$ in the upper left, $x,y$ adjacent to it, and $xy$ in the lower right.
Remember you only get one edge per generator. You do not get edges for anything you can get to eventually, otherwise your graph would always be a complete graph by the definition of a generating set.


text:
\section{Introduction}
\label{Sec: Introduction}

Fractional models construct a tractable mathematical framework to describe and predict the behavior of multi-scales multi-physics complex phenomena. Particularly, fractional differential equations, as a well-structured generalization of their integer order counterparts, provide a rigorous mathematical tool to develop models, describing anomalous behavior in complex physical systems \cite{zhang2017review,jaishankar2014,jha2003evidence, Castillo2004Plasma,jaishankar2013,naghibolhosseini2015estimation,naghibolhosseini2017fractional,meer01, meral2010fractional}, where the anomalies manifest in heavy tailed distribution of corresponding underlying stochastic processes, moreover, exhibiting sharp peaks, intermittency, and asymmetry in the underlying distributions. Significant approximations as inherent part of assumptions upon which the model is built, lack of information about true values of parameters (incomplete data), and random nature of quantities being modeled pervade uncertainty in the corresponding mathematical formulations \cite{cullen1999probabilistic,roy2011comprehensive}. In this work, we develop an uncertainty quantification (UQ) framework in the context of stochastic fractional partial differential equations (SFPDEs), in which we characterize different sources of uncertainties and further propagate the associated randomness to the system response quantity of interest (QoI).

\vspace{0.1 in}

\noindent\textbf{Types and Sources of Uncertainty}.

The model uncertainties are in general being classified as aleatory, epistemic, and mixed, according to their fundamental essence. They can also be broadly characterized as occurring in model inputs, numerical approximations, and model form. Model inputs encompass all model parameters coming from geometry, constitutive laws, and fields equation, while also pertaining surrounding interactions, such as boundary conditions and random forcing sources (noise). Numerical approximations, which are an essence of differential equations since they generally do not lend themselves to analytical solutions, introduce uncertainty by imposing different sources of discretization error, iterative convergence error, and round off error. The fractional derivatives introduce derivative orders, namely fractional indices, as new set of model parameters in addition to model coefficients. These new parameters are strongly tied to the distribution of underlying stochastic process and their statistics are estimated from experimental observations in practice, see e.g.  \cite{benson2000application,baeumer2001subordinated}.

\vspace{0.1 in}

\noindent\textbf{Uncertainty Framework}.

Conventional approaches in parametric UQ of differential equations is based around Monte Carlo sampling (MCS) \cite{carlo1996concepts}, which performs ensemble of forward calculations to map the uncertain input space to the uncertain output space.
This method enjoys from being embarrassingly parallelizable and can be implement quite readily on high dimensional random spaces. However, the key issue is the slow rate of convergence $\sim 1/\sqrt{N}$ with $N$ number of realization, which consequently imposes exhaustively so many operations of forward solver, makes it not practical for expensive simulations. Other methods such as sequential MCS \cite{del2006sequential} and multilevel sequential MSC \cite{beskos2017multilevel} are also developed and recently used in \cite{jasra2016forward} to improve the parametric uncertainty assessment in elliptic nonlocal equations. An alternative to expensive MCS is to build surrogate models. An extensive comparison of two widely used ones, namely polynomial chaos and Gaussian process, are provided in a recent work \cite{owen2017comparison}. Polynomial chaos, in which the output of stochastic model is represented as a series expansion of input parameters, was initially applied in \cite{ghanem2003stochastic}, and later extended and used in \cite{xiu2002modeling,xiu2002wiener,knio2006uncertainty}. It is also generalized and used in constructing stochastic Galerkin methods \cite{babuska2004galerkin,babuvska2005solving,le2004uncertainty,le2004multi} for problems with higher-dimensional random spaces. Other non-sampling numerical methods, including but not limited to perturbation method \cite{schuss1980singular,babuvska2002solving,todor2005sparse,winter2002groundwater} and moment equation method \cite{liu1986probabilistic,liu1986random} are also developed, however their applications are restricted to stochastic systems with relatively low-dimensional random space. These so-called ``\textit{intrusive}" approaches typically do not treat the forward solver as a black-box, rather require some knowledge and reformulation of the governing equations and thus, may not be practical in many problems with complex codes.

A wide range of ``\textit{non-intrusive}" techniques mostly stretch over sampling, quadrature, and regression, see \cite{owen2017comparison} and references therein. More recently, high-order probabilistic collocation methods (PCM), employing the idea of interpolation/collocation in the random spaces, are developed in \cite{xiu2005high,babuvska2007stochastic,nobile2008sparse}. These methods limit the sample points to an efficient subset of random space, while adequately sampling the necessary range. The excellence in use of PCM is twofold; it has the benefit of easily sampling at nodal points that naturally leads to independent realizations of the problem as in MCS, and the advantage of fast convergence rate. The challenging post processing of solution statistics, which can essentially be described as a high-dimensional integration problem, can also be resolved by adopting sparse grid generators, such as Smolyak algorithm \cite{nobile2008sparse,Smolyak63}. The use of sparse grids, as opposed to full tensor product construction from one-dimensional quadrature rules, will effectively reduce the number of sampling, while preserving a fast convergence rate to high level of accuracy.

\vspace{0.1 in}

\noindent\textbf{Forward Solver}.

A core task in computational forward UQ is to form an efficient numerical method, which for each realizations of random variables can accurately solves and simulates the deterministic counterpart of stochastic model in the physical domain. Such numerical method is usually called ``\textit{forward solver}" or ``\textit{simulator}". In the case of FPDEs, the excessive cost of numerical approximations due to the inherent nonlocal nature of fractional differential equations additionally become more challenging when generally most of uncertainty propagation techniques instruct operations of forward solver many times. This requires implementation of more efficient numerical schemes, which can manage increasing computational costs while maintaining sufficiently low error level in mitigating the corresponding uncertainties. In addition to numerous finite difference methods for solving FPDEs \cite{Gorenflo2002, Sun2006,  Lin2007, wang2010direct, wang2011fast, Cao2013, zeng2015numerical, zayernouri2016_JCP_Frac_AB_AM}, recent works have elaborated efficient spectral schemes, for discretizing FPDEs in physical domain, see e.g., \cite{Rawashdeh2006, Lin2007, Khader2011, Khader2012, Li2009, Li2010, chen2014generalized, wang2015high, bhrawy2015spectral}. More recently, Zayernouri et al. \cite{Zayernouri2013, zayernouri2015tempered} developed two new spectral theories on fractional and tempered fractional Sturm-Liouville problems, and introduced explicit corresponding eigenfunctions, namely \textit{Jacobi poly-fractonomials} of first and second kind. These eignefunctions are comprised of smooth and fractional parts, where the latter can be tunned to capture singularities of true solution. They are successfully employed in constructing discrete solution/test function spaces and developing a series of high-order and efficient Petrov-Galerkin spectral methods, see \cite{lischke2017petrov, suzuki2016fractional,samiee2016,samiee2017fast,samiee2017unified,samiee2018petrov,kharazmi2017petrov,kharazmi2017sem,kharazmi2018fractional}.

\vspace{0.1 in}

The main focus of this work is to develop an operator-based computational forward UQ framework in the context of stochastic fractional partial differential equation. Assuming that the mathematical model under consideration is well-posed and accounts in principle for all features of underlying phenomena, we identify three main sources of uncertainty, \textit{i}) parametric uncertainty, including fractional indices as new set of random parameters appeared in the operator, \textit{ii}) additive noises, which incorporates all intrinsic/extrinsic unknown forcing sources as lumped random inputs, and \textit{iii}) numerical approximations. Computational challenges arise when the admissible space of random inputs is infinite-dimensional, e.g. problems subject to additive noise \cite{rizzi2012uncertainty}, and thus, the framework involves uncertainty parametrization via a finite number of random space basis. Unlike the classical approach in modeling random inputs, which considers idealized uncorrelated processes (white noises), we model the random inputs as more/fully correlated random processes (colored noises), and parametrize them via Karhunen-Lo\`{e}ve (KL) expansion by assuming finite-dimensional noise assumption. This yields the problem in finite dimensional random space. We then, propagate the parametric uncertainties into the system response by applying PCM. We obtain the corresponding deterministic FPDE for each realization of random variables, using the Smolyak sparse grid generators for low to moderately high dimensions. In order to formulate the forward solver, we follow \cite{samiee2016} and develop a high-order Petrov-Galerkin (PG) spectral method to solve for each realization of SFPDE in the physical domain, employing Jacobi poly-fractonomials in addition to Legendre polynomials as temporal and spatial basis/test functions, respectively. The smart choice of coefficients in construction of spatial basis/test functions yields symmetric properties in the resulting mass/stiffness matrix, which is then exploited to formulate an efficient fast solver. We also show that for each realization of random variables, the deterministic problem is mathematically well-posed and the proposed forward solver is stable. By adopting sufficient number of basis in the physical domain, we eliminate the epistemic uncertainty associated with numerical approximation and isolate the impact of parametric uncertainty on system response QoI.

The organization of the paper is as follows. We recall some preliminaries on fractional calculus in section \ref{Sec: Fractional Calculus}. Then, we formulate the stochastic system in section \ref{Sec: Uncertainty framework}, and parametrize the random inputs. We also develop the stochastic sampling, namely PCM and MCS for our stochastic problem. We further construct the deterministic solver in section \ref{Sec: Forward Solver}, and provide the numerical results in section \ref{Sec: numerical results}. We end the paper with a conclusion and summary.

config: arxiv


TITLE: Percentage question based on voting
QUESTION [2 upvotes]: Please help me solve this problem. Thank you.

REPLY [1 votes]: Assume that Mondale has 60 supporters and Bush 40. Then 6 supporters changed to Bush and 4 to Mondale. So, now Mondale has 60+4-6=58 supporters and Bush has 40+6-4=42 supporters. You need to find the number of supporters $x$ expressed as a percentage of the total number of voters (100), so that 58 becomes 50 (at par). So, $8$ in $58$ is equal to $x$ in 100. Solve for $x$:
$$x=\frac{8}{58}\cdot100\approx14\%$$