Singular Value Decay

Approximting Eigenvalues

Definition (Moments $\mu_k$). Let $\Gamma$ be a positive oriented closed Jordan curve in the complex plane. Let $\lambda_1,\ldots,\lambda_m$ be distinct eigenvalues that lie in the interior of $\Gamma$. Then $$ \mu_k = \frac{1}{2\pi i}\int_\Gamma (z-\gamma)^k f(z)\ dz = \sum_{j=1}^ma_j(\lambda_j-\gamma)^k$$ where $\gamma$ is located inside $\Gamma$.

Define the Hankel matrices $$ H = \left[\begin{array}{cccc} \mu_0 & \mu_1 & & \mu_{m-1}\\ \mu_1 & \mu_2 & & \mu_{m}\\ & & \ddots& \\ \mu_{m-1} & \mu_{m} & & \mu_{2m-2}\\ \end{array}\right] \quad\mbox{ and }\quad H_s = \left[\begin{array}{cccc} \mu_1 & \mu_3 & & \mu_{m}\\ \mu_2 & \mu_3 & & \mu_{m+1}\\ & & \ddots& \\ \mu_{m} & \mu_{m+1} & & \mu_{2m-1}\\ \end{array}\right]$$ Let $V_m$ be the Vandermonde matrix $$V_m= \left[\begin{array}{cccc} 1 & 1 & & 1\\ \lambda_1-\gamma & \lambda_2-\gamma & & \lambda_m-\gamma\\ & & \ddots& \\ (\lambda_1-\gamma)^{m-1} & (\lambda_2-\gamma)^{m-1} & & (\lambda_m-\gamma)^{m-1}\\ \end{array}\right]$$ Assuem $a_i\neq0$ for $i = 1,2,\ldots,m$ and define $$D_m = \texttt{diag}(a_1,\ldots,a_m)\quad\mbox{and}\quad \Lambda_m = \texttt{diag}(\lambda_1-\gamma,\ldots,\lambda_m-\gamma) $$ Therefore, we can verify $$ H = V_m D_m V_m^T \quad\mbox{and}\quad H_s = V_m D_m \Lambda_m V_m^T$$ Since $$ H = V_m D_m V_m^T \quad\mbox{and}\quad H_s = V_m D_m \Lambda_m V_m^T$$ Therefore $$\begin{array}{rcl} H_s - \lambda H &=& V_m D_m\Lambda_m V_m^T - \lambda V_m D_m V_m^T \\ \\ & =& V_m D_m (\Lambda_m - \lambda I ) V_m^T \end{array}$$ Since $V_m, D_m, V_m$ are nonsingular, then the eigenvalues of the pencil $H_s - \lambda H$ are given by $\lambda_i-\gamma$ for $i =1,2,\ldots,m$.\\ Therefore, we can obtain the eigenvalues $\lambda_1,\ldots,\lambda_m$ by solving the generalized eigenvalue problem $H_s x = \lambda H x$. If we find an appropriate $\Gamma$ that includes a small number of eigenvalues, then the size of the derived eigenproblem is small.

This ideas are implemented as the Sakurai-Siguira algorithm whcih lead to good results in general. alt text

Algebraic Multigrid

As finite element discretizations ever grow in size to address real word problems, there is an increasing need\n for fast algorithms. Nowadays there are many GPU/CPU parallel approaches to solve such problems. Multigrid methods can be used to solve such large-scale problems, or even better they can be used to precondition the conjugate gradient method, yielding better results in general.

Capabilities of multigrid algorithms relies on on the effectiveness of the inter-grid transfer operators. Multigrid methods exploit discretizations of a given problem using different mesh sizes, seeking better convergence than iterative methods. These techniques are founded in the principle of divide and conquer.

alt text

Iterative methods do a good job of damping high-frequency modes, but low-frequency modes remains almost the same, converging very slow. But there low-frequency modes can be represented well on coarse meshes, which are low in dimension and thus easier to solve. Therefore more error components can be eliminated by transferring the problem from a fine mesh to a coarse mesh. Obviously, this process can be repeated using a hierarchy of meshes or grids.

Aggregation Graph Partitioning

Given a graph $G = (V;E)$, where $V$ is the set of vertices (or nodes), and $E$ is the set of edges. An edge is a pair of two distinct vertices. Then, MIS($k$) is defined as follow.

Definition (MIS($k$)). Given a graph $G = (V;E)$, let $V_{root} \subseteq V$ be a set of root vertices, and denote by $dG(\cdot;\cdot)$ the length of the shortest path between two vertices in the graph. Then $V_{root}$ is a maximal independent set of distance $k$, or MIS($k$), if the following both hold.

  1. Independence: Given any two vertices $u,v \in V_{root}$, then $dG(u; v) > k$.
  2. Maximality: There exists no $u \in V\setminus V_{root}$ such that $dG(u; v) > k$ for all $v \in V_{root}$.

Given a valid MIS(2), the aggregates are constructed as follows: assume the MIS(2) is specified by an array of values $\{0;1\}$, where a value of 1 means the corresponding node is a member of the MIS(2) and a value 0 otherwise. Since the $k$-th node in the MIS(2) serves as the root of the $k$-th aggregate, the only remaining task is to propagate the root indices outward. For each node $j$ not in the MIS(2), we evaluate the first ring of nodes strongly connected to $j$. If there is a root node among these nodes, we will assign the node $j$ to the most strongly connected root node. Then we evaluate the second ring of nodes that are strongly connected to $j$; then the node $j$ is assigned to the aggregate whose root node is more strongly connected to $j$ within the first and second ring of strongly connected nodes.


Wavelet Decomposition/Reconstruction

A multiresolution analysis (MRA) or multiscale approximation (MSA) is the design method of most of the practically relevant discrete wavelet transforms (DWT) and the justification for the algorithm of the fast wavelet transform (FWT). It was introduced in this context in the 80s by Mallat and Meyer.

Definition (MRA). An MRA is a sequence $\{V_j\}_{j\in\mathbb{Z}}$ of subspaces of $L^2(\mathbb{R})$ such that

  1. $V_j \subset V_{j+1}$ for all $j\in\mathbb{Z}$.
  2. $\displaystyle\overline{\bigcup_{j\in\mathbb{Z}}V_j } = L^2(\mathbb{R})$.
  3. $\displaystyle\bigcap_{j\in\mathbb{Z}} V_j = \{0\}$.
  4. $f \in V_j$ if and only if the function $x \mapsto f(2^{-j}x) \in V_0$.
  5. $f\in V_0$ if and only if the funciton $x\mapsto f(x-m) \in V_0$ for all $m\in\mathbb{Z}$.
  6. There exists a function $\phi\in V_0$, called scaling function, such that $\{ \phi(t-m) \}_{m\in\mathbb{Z}}$ is an orthonormal basis of $V_0$.

Given an MRA $\{V_j\}_{j\in \mathbb{Z}}$ each subspace $V_j$ can be expresed as $V_j = V_{j-1}\oplus W_{j-1}$. Let $P_j$ and $Q_j$ the orthogonal projectors onto $V_j$ and $W_j$ respectively. Thus each element $f\in V_j$ can be projected into a lower resolution space $V_{j-1}$ by $P_jf$, this proces could be done iteratively leading to the decomposition algorithm.

Using a threshold filtering over the results we could get different results applying the reconstruction algorithm.