{"id":723,"date":"2024-04-01T21:31:42","date_gmt":"2024-04-01T20:31:42","guid":{"rendered":"https:\/\/sinatootoonian.com\/?p=723"},"modified":"2025-12-27T16:05:07","modified_gmt":"2025-12-27T16:05:07","slug":"notes-on-multiresolution-matrix-factorization","status":"publish","type":"post","link":"https:\/\/sinatootoonian.com\/index.php\/2024\/04\/01\/notes-on-multiresolution-matrix-factorization\/","title":{"rendered":"Notes on Multiresolution Matrix Factorization"},"content":{"rendered":"\n

These are my notes from early January on Kondor et al.’s Multiresolution Matrix Factorization<\/a> from 2014. This was a conference paper and the exposition was a bit terse in places, so below I try to fill in some of the details I thought were either missing or confusing. <\/p>\n\n\n\n

Motivating MMF<\/h3>\n\n\n\n

We will be interested in factorizing $n \\times n$ symmetric matrices $A$. MMF can be motivated through PCA, which gives the eigendecomposition $A = U^T D U,$ where $U$ is orthonormal and $D$ is diagonal. We can think of this as a compression if $U$ is known in advance. We can then encode $A$ (and other such matrices) by their projections onto $U$. For $A$, this is simply $$ U A U^T = D.$$ But we can imagine dealing with families of matrices and storing only their projections onto a fixed $U$. If $U$ is adapted to the family then these projections will be sparse, just like the diagonal matrix $D$ was (at most $n$ of the $n^2$ elements non-zero), and we will have achieved compression. <\/p>\n\n\n\n

The problem with PCA is that the basis elements $U$ are non-local. That is, a given basis element will have lots of non-zero entries, implying that the corresponding projection requires information from many of the $n$ locations. This can make the projections hard to interpret.<\/p>\n\n\n\n

MMF takes ideas from multiresolution analysis of functions on the real line and applies them to the discrete setting of $n$-dimensional vectors. It performs the decomposition of $A$ sequentially using a series of sparse orthonormal matrices. So our decomposition of $A$ will look like
$$ U_L U_{L-1} \\dots U_{1} A U_1^T \\dots U_{L-1}^T U_{L} = H,$$ where $H$ will also be diagonal except possibly on some central block. That is, we decompose the global rotation $U$ that a procedure like PCA might produce into a product of a sequence of local rotations $U_1, U_2, \\dots U_L$.<\/p>\n\n\n\n

The MMF Objective<\/h3>\n\n\n\n

To find the decomposition above we solve the optimization problem \\begin{align} \\min_{ \\substack{[n]\\supseteq S_1\\dots S_L \\\\ U_1,\\dots,U_L \\in O \\\\ H \\in H^n_{S_L}}} \\|U_L \\dots U_1 A U_1^T \\dots U_L^T – H \\|_F^2. \\label{opt}\\tag{1}\\end{align}<\/p>\n\n\n\n

Here the $[n]=S_0 \\supseteq S_1 \\dots S_{L-1} \\supseteq S_L$ are nested subsets of the integers $1 \\dots n$ that indicate the coordinates that can still be involved in rotations. All remaining coordinates are held fixed. $H_{S_L}^n$ is the set of matrices that are non-zero on the $S_L \\times S_L$ block and also on the diagonal. <\/p>\n\n\n\n

The first step is to eliminate the last constraint in the minimization ($H \\in H_{S_L}^n$) by moving it to the norm itself. We do this by defining a ‘residual’ norm that only counts elements outside of the block and diagonal: \\begin{align} \\|X\\|^2_\\text{resi} = \\sum_{\\substack{(i,j) \\not \\in S_L \\times S_L, \\\\ i \\neq j }} X_{ij}^2. \\label{resi} \\tag{2}\\end{align} So our problem becomes \\begin{align} \\min_{\\substack{[n]\\supseteq S_1\\dots S_L \\\\ U_1,\\dots,U_L \\in O}} \\|U_L \\dots U_1 A U_1^T \\dots U_L^T \\|_\\text{resi}^2. \\label{resi2} \\tag{3}\\end{align} In other words, we’re looking for a sequence of rotations that diagonalize $A$, except possibly on the $S_L \\times S_L$ block. <\/p>\n\n\n\n

Multiresolution analysis of $\\mathbb R^n$.<\/h3>\n\n\n\n

We now need to be careful about what the $U_i$ are. We first express $A$ as $$ A = \\sum \\la_i \\xx_i \\xx_i^T.$$ Each of the $\\xx_i \\in \\RR^n = V_0$. The idea is now to split the base space $V_0$ into a direct sum of a ‘smooth’ subspace $V_1$ and a ‘detail’ subspace $W_1$: $$ V_0 = V_1 \\oplus W_1.$$ <\/p>\n\n\n\n

Let a basis for $V_1$ be the rows of $\\Phi_1$ and a basis for $W_1$ the rows of $\\Psi_1$. We can express these basis elements in terms of the basis $\\Phi_0$ for $V_0$ using the matrix $U_1$: $$ \\begin{bmatrix} \\Phi_1 \\\\ \\Psi_1 \\end{bmatrix} = U_1 \\Phi_0.\\label{basis0}\\tag{4}$$ Since $\\dim(\\Phi_1) = \\dim(V_1)$ and $\\dim(\\Psi_1) = \\dim(W_1)$ and $\\dim(V_1) + \\dim(W_1) = \\dim(V_0)$, $U_1$ has dimensions $\\dim(V_0) \\times \\dim(V_0) = n \\times n$. We additionally require that $U_1$ is orthonormal, i.e. $U_1^T U_1 = I$. <\/p>\n\n\n\n

We can now express $A$ in the $\\Phi_1 \\cup \\Psi_1$ basis by first using $\\Eqn{basis0}$ to express $\\Phi_0$ in this basis: \\begin{align} U_1^T \\begin{bmatrix} \\Phi_1 \\\\ \\Psi_1 \\end{bmatrix} = \\Phi_0 \\implies \\Phi_0^T = \\begin{bmatrix} \\Phi_1^T | \\Psi_1^T \\end{bmatrix} U_1.\\end{align} Then $$ \\xx_i = \\Phi_0^T \\xx_i = \\begin{bmatrix} \\Phi_1^T | \\Psi_1^T \\end{bmatrix} U_1 \\xx_i.$$ The first $\\dim(V_1)$ rows of $U_1$, whose projections on $\\xx_i$ are projected out with $\\Phi_1^T,$ are then the scaling functions, and the remaining rows are the wavelets. <\/p>\n\n\n\n

\nSo the coefficients expressed in the new basis are $U_1 \\xx_i$. Therefore,\n$$ A_1 = \\sum \\la_i U_1 \\xx_i (U_1 \\xx_i)^T = U_1 \\left( \\sum \\la_i \\xx_i \\xx_i^T \\right) U_1^T = U_1 A U_1^T.$$\n<\/p>\n\n\n\n

To continue with the multiresolution analysis, we split $V_1 = V_2 \\oplus W_2$ and define $U_2$ as the matrix that expresses $\\Phi_2$ and $\\Psi_2$ in terms of $\\Phi_1$: $$ \\begin{bmatrix} \\Phi_2 \\\\ \\Psi_2 \\end{bmatrix} = U_2 \\Phi_1.$$ The dimensions of $U_2$ are $\\dim(V_1) \\times \\dim(V_1)$, and as before we also require that $U_2$ is orthogonal, i.e. $U_2^T U_2 = I$. To make it compatible with $A_1$, we extend it with the identity matrix, so that it leaves the $W_1$ subspace unchanged: $$ U_2 \\leftarrow U_2 \\oplus I_{\\dim(W_1)}.$$ Then, following the same reasoning as before, we can express $A$ in the $\\Phi_2 \\cup \\Psi_2 \\cup \\Psi_1$ basis as $$ A_2 = U_2 A_1 U_2^T = U_2 U_1 A U_1^T U_2^T.$$ The first $\\dim(V_2)$ rows of $U_2 U_1$ are the scaling functions, the next $\\dim(W_2)$ rows are the wavelets from $W_2$, and the remaining rows are the wavelets from $W_1.$ Continuing in this way for $L$ steps, we get the optimization problem $\\Eqn{opt}.$ <\/p>\n\n\n\n

Recursive decomposition<\/h3>\n\n\n\n

To solve our optimization problem $\\Eqn{opt}$ we work with its residual-norm formulation $\\Eqn{resi2}$ and expand it recursively. First, let’s apply a recursive decomposition to the ordinary Frobenius norm. Since $S_0 = [n]$, we start with $$ \\|H\\|^2 = \\|H_{S_0, S_0}\\|^2$$. This is the sum of squares of all the elements of $H$. The next set of indices, $S_1$ are a proper subset of $S_0$, so we can decompose our sum as
$$ \\|H_{S_0, S_0}\\|^2 = \\|H_{S_1, S_1}\\|^2 + \\|H_{S_1, S_0\/S_1}\\|^2 + \\|H_{S_0\/S_1,S_1}\\|^2 + \\|H_{S_0\/S_1,S_0\/S_1}\\|^2.$$ Letting $J_i = S_{i-1}\/S_i$ indicate the indices $S_{i-1}$ contains that $S_i$ doesn’t, we can write the last expression as $$ \\|H_{S_0, S_0}\\|^2 = \\|H_{S_1, S_1}\\|^2 + \\|H_{S_1, J_1}\\|^2 + \\|H_{J_1,S_1}\\|^2 + \\|H_{J_1,J_1}\\|^2.$$ We can pack all the terms that index $J_1$ into $\\mathcal E_1$ and get $$ \\|H_{S_0, S_0}\\|^2 = \\|H_{S_1, S_1}\\|^2 + \\mathcal{E}_1.$$<\/p>\n\n\n\n

We can similarily decompose $\\|H_{S_1,S_1}\\|^2$, to get \\begin{align}\\|H_{S_0, S_0}\\|^2 &= (\\|H_{S_2, S_2}\\|^2 + \\|H_{S_2, J_2}\\|^2 + \\|H_{J_2,S_2}\\|^2 + \\|H_{J_2,J_2}\\|^2) + \\mathcal{E}_1\\\\ &= \\|H_{S_2, S_2}\\|^2 + \\mathcal{E}_2 + \\mathcal{E}_1.\\end{align} Continuing in this way, we get $$\\|H\\|^2 = \\|H_{S_0, S_0}\\|^2 = \\|H_{S_L, S_L}\\|^2 + \\sum_{\\ell=1}^L \\mathcal{E}_\\ell.$$<\/p>\n\n\n\n

Therefore, if we redefine $\\mathcal{E}_i$ to not include the diagonal terms, $$\\mathcal{E}_\\ell = 2 \\|H_{J_\\ell, S_{\\ell}}\\|_F^2 + \\|H_{\\substack{J_\\ell, J_\\ell\\\\ \\text{off diag}}}\\|^2, \\label{Eell} \\tag{5}$$ we arrive at $$\\|U_L \\dots U_1 A U_1^T \\dots U_L^T \\|_\\text{resi}^2 = \\|H\\|_\\text{resi}^2 = \\sum_{\\ell=1}^L \\mathcal{E}_\\ell.$$ <\/p>\n\n\n\n

We’d now like to express $\\Eqn{Eell}$ succinctly in terms of $A$ and the $U$’s. To do so, notice that $U_{\\ell+1}$ is the identity on its lower $S_0\/S_{\\ell}$ block. So, with some abuse of notation, \\begin{align} A_{\\ell+1} = U_{\\ell+1} A_{\\ell} U_{\\ell+1}^T &= \\begin{bmatrix} U_{\\ell+1} & 0 \\\\ 0 & I_{S_0\/S_{\\ell}} \\end{bmatrix} \\begin{bmatrix} A_{11} & A_{12} \\\\ A_{21} & A_{22} \\end{bmatrix} \\begin{bmatrix} U_{\\ell+1}^T & 0 \\\\ 0 & I_{S_0\/S_{\\ell}} \\end{bmatrix}\\\\ &= \\begin{bmatrix} U_{\\ell+1} A_{11} & U_{\\ell+1} A_{12} \\\\ A_{21} & A_{22} \\end{bmatrix} \\begin{bmatrix} U_{\\ell+1}^T & 0 \\\\ 0 & I_{S_0\/S_{\\ell}} \\end{bmatrix}\\\\ &= \\begin{bmatrix} U_{\\ell+1} A_{11} U_{\\ell+1}^T & U_{\\ell+1} A_{12} \\\\ A_{21} U_{\\ell+1}^T & A_{22} \\end{bmatrix}. \\end{align} Then $$ \\|[A_\\ell]_{J_\\ell, S_\\ell}\\|_F^2 = \\|[A_\\ell]_{S_\\ell, J_\\ell}\\|_F^2 = \\|[A_{12}]_{:,J_\\ell}\\|_F^2,$$ where we’ve abused notation on the final set of indices and used the absolute indices $J_\\ell$ instead of those relative to $A_{12}$. Next, $$ \\|[A_{\\ell+1}]_{J_\\ell, S_\\ell}\\|_F^2 = \\|[A_{\\ell+1}]_{S_\\ell, J_\\ell}\\|_F^2 = \\|[U_{\\ell+1} A_{12}]_{:,J_\\ell}\\|_F^2 = \\|[A_{12}]_{:,J_\\ell}\\|_F^2 = \\|[A_\\ell]_{S_\\ell, J_\\ell}\\|_F^2,$$ where the third equality is because Frobenius norm is invariant to rotation. But this will then hold for all $\\ell + k$, because $U_{\\ell+k}$ is a rotation on a subspace of that on which $U_{\\ell+1}$ acts. Therefore $$ \\|H_{J_\\ell, S_\\ell}\\|_F^2 = \\|[A_L]_{J_\\ell, S_\\ell}\\|_F^2 = \\|[A_\\ell]_{J_\\ell, S_\\ell}\\|_F^2.$$ <\/p>\n\n\n\n

The only term that remains in $\\Eqn{Eell}$ is $\\|H_{J_\\ell, J_\\ell}\\|_F^2$ (without the diagonal terms) – but these will be part of $A_{22}$, which is left unchanged in $A_{\\ell+1}$. Therefore, we can write $$ \\mathcal{E}_\\ell = 2 \\|[A_\\ell]_{J_\\ell, S_\\ell}\\|_F^2 + \\|[A_\\ell]_{\\substack{J_\\ell, J_\\ell \\\\ \\text{off diag}}}\\|_F^2.$$ <\/p>\n\n\n\n

\nNotice that $S_\\ell$ is potentially a large set, while the rotations might be operating on small subsets of this. The rotation $U_\\ell$ can operate on any of $S_{\\ell-1}$. Let’s suppose that: \n<\/p>\n\n\n\n