Skeleton decomposition

Digital images, integral and differential operators are usually stored in matrix form. These matrices are often very large systems and require enormous memory and processing power. In numerical linear algebra, mathematicians are interested in developing algorithms that can efficiently reduce storage without losing too much accuracy. This field is known as low rank approximation.

Many matrices that arise from partial differential equations or integral equations are low rank systems because of translational invariance. The most famous and probably most well known low rank approximation algorithm is the singular value decomposition method. Given a matrix A\in R^{m\times n} , we can decompose A into the product of three matrices in the form of

A=USV^{T},

where U\in m\times r , S\in r\times r , and V\in n\times r . In this setting, S is a diagonal matrix and its entries are called the singular values. By convention, the singular values are listed in decreasing order along the diagonal of S. In many systems, the relative magnitude of singular values decays very quickly:

\frac{\sigma_{i}}{\sigma_{1}}\rightarrow0

for large i . For this reason, we can approximate the matrix A by only using the first few singular values. We can write

A\simeq R=\hat{U}\hat{S}\hat{V}^{T} ,

where \hat{U}\in m\times r^{*} ,\hat{S}\in r^{*}\times r^{*} ,\mbox{and }\hat{V}\in n\times r^{*},r^{*}\ll r . One can prove that singular value decomposition yields the best possible low rank approximation in the Frobenius norm, however this method is computational costly and sometimes numerically unstable. Instead, we introduce a much more efficient method: skeleton decomposition.

The idea is fairly simple, we can rewrite A as

A\backsimeq C\hat{A}^{-1}R ,

where C is the column restriction of A , R is the row restriction of A, and \hat{A} is the intersection of row and column restriction. In another word, if A has rank r, then we pick r columns from A and form a submatrix, and then we pick r rows from A and form a submatrix. Finally, the intersection of C and R gives \hat{A} . The figure below illustrates this process.

Image

Skeleton decomposition: A\backsimeq C\hat{A}^{-1}R ,

This method does not need much computation, since we are simply picking rows and columns! The question is how do we pick the best rows and the best columns? This is actually a quiet interesting optimization problem. It turns out that we want to pick the rows and columns that maximize the spanned volume. In another words, we maximize |\det(\hat{A})| . Now this question looks like the maximum ellipsoid problem we have on problem set 3. There are multiple iterative algorithms that search for the optimal set of rows and columns. Look into adaptive/ incomplete cross approximation methods for details.

Reference

http://math.mit.edu/icg/papers/sublinear-skeleton.pdf

Written by Yufeng (Kevin) Chen

Leave a comment