Learning to Solve Multiresolution Matrix Factorization by Manifold Optimization and Evolutionary Metaheuristics

We frame the task of identifying the optimal nested sequence of sets ${\mathbb{S}}_{L}\subseteq\cdots\subseteq{\mathbb{S}}_{1}\subseteq{\mathbb{S}}_{0}=[n]$ as learning a set of wavelet indices to solve the MMF optimization in (2). This involves two primary components for index selection at each resolution level $\ell\in\{1,\ldots,L\}$:

• •

  Select the set of indices ${\mathbb{T}}_{\ell}\subset{\mathbb{S}}_{\ell-1}$ representing the rows/columns to be wavelets at this level, which are then eliminated by defining ${\mathbb{S}}_{\ell}={\mathbb{S}}_{\ell-1}\setminus{\mathbb{T}}_{\ell}$. To simplify computation, we assume that each resolution level selects only one row/column as the wavelet, such that $\lvert{\mathbb{T}}_{\ell}\rvert=1$. Consequently, the cardinality of ${\mathbb{S}}_{\ell}$ decreases by 1 at each level, giving $d_{\ell}=n-\ell$. The core block size of ${\bm{H}}$ becomes $(n-L)\times(n-L)$, corresponding to exactly $n-L$ active rows/columns at the end.

• •

  Select $k-1$ indices ${\mathbb{I}}_{\ell}=\{i_{1},\ldots,i_{k-1}\}\subset{\mathbb{S}}_{\ell-1}$ to construct the corresponding rotation matrix ${\bm{U}}_{\ell}$ (see Section  4).

A small example to illustrate this index selection is given in Fig. 1. The two algorithms use the Frobenius error from the Stiefel manifold optimization as the fitness function. In the second step of selecting $\ignorespaces k-1$ indices, we identify the $\ignorespaces k-1$ indices whose rows are closest to the wavelet row in Euclidean distance (see 11.3 for the inspiration of this heuristic), reducing this problem to finding an ordered set of $\ignorespaces L$ wavelet indices.

In both metaheuristics, the candidate solution is an ordered set of $\ignorespaces L$ indices chosen as wavelet indices. The fitness function employed is the initial cost from the Stiefel manifold optimization algorithm, without training iterations, to estimate solution quality. This approach is designed to minimize the computational time when running the metaheuristics, as the optimization phase of the Stiefel manifold algorithm is only executed after identifying the best solution through the metaheuristics. Thus, the initial cost from the Stiefel manifold optimization serves as a cost-effective estimate of solution quality.

We employ a metaheuristics-based approach grounded in evolutionary algorithms to solve the optimization problem of finding the optimal nested sequence for Multiresolution Matrix Factorization (MMF).

Evolutionary algorithms [33], inspired by the process of natural selection and genetics, are particularly effective for complex optimization tasks. Our method iteratively improves a population of candidate solutions by applying operations such as selection, crossover, and mutation. The selection process identifies the most promising candidates based on a fitness function, while crossover and mutation introduce genetic diversity, enabling the exploration of the solution space.

The mutation part of the algorithm consists of two independent mutation operators. The first mutation operator is randomly swapping two indices of the candidate solution. The second mutation operator is replacing a random value in the ordered set with another value which is not already in the solution. A solution cannot have duplicated values.

The crossover operator used in our approach is a random one-point crossover. However, a challenge arises because this crossover can create invalid offspring with duplicated elements if both parents have common elements. To address this issue, we implement a strategy to separate the values common to both parents from the other values.

Specifically, the values that are common to both parents will be preserved and not subjected to crossover. Next, we separate the remaining values (those not common to both parents) into two new sets of genes. These new sets will contain only unique values, ensuring that no duplicates are present. A normal one-point crossover can be perfomed on these new gene sets, creating two new sets of genes without any common values. After the crossover, reinsert the common values back into the respective offspring. This ensures that the offspring are valid and maintain the necessary elements from both parents without any duplication. Fig. 2 details an example using this crossover operator.

By using this method, we ensure that the resulting offspring are valid and retain genetic diversity from both parents while avoiding any duplicate values.

Directed evolution, a laboratory methodology wherein biological entities possessing desired characteristics are generated through iterative cycles of genetic diversification and library screening or selection, has emerged as a highly valuable and extensively utilized instrument in both fundamental and practical realms of biological research [34, 35, 36].

This directed evolution algorithm uses the same mutation operators as the evolutionary algorithm. Mutation is performed on every member of the parent population. The parent population is part of the next generation.

Figure 3 illustrates the convergence behavior of two metaheuristics applied to the Karate Club matrix, contrasted with various random and heuristic baselines. Notably, after 100 generations, the Evolutionary Algorithm (EA) surpasses all other baselines in performance. Additionally, Directed Evolution (DE) demonstrates effective MMF approximation. However, it is worth noting that EA requires significantly more time to reach convergence over 100 generations. This disparity in time consumption is attributed to the larger population size initialized in EA compared to DE. Other heuristic baselines, although quicker, fail to provide comparably accurate approximations as the two metaheuristics.

Figure 4 depicts the wavelet bases at different levels of resolution. The low index wavelets are localized since it teases out distinct local structures of the matrix, separating it into rough and smoother components. The higher wavelet bases are left only with the smooth part of the matrix. This observation aligns with the interpretation provided for MMF in 11.2.

The eigendecomposition of the normalized graph Laplacian operator $\tilde{{\bm{L}}}={\bm{U}}^{T}{\bm{H}}{\bm{U}}$ can be used as the basis of a graph Fourier transform. [37] defines graph Fourier transform (GFT) on a graph $\mathcal{G}=(V,E)$ of a graph signal ${\bm{f}}\in\mathbb{R}^{n}$ (that is understood as a function $f:V\rightarrow\mathbb{R}$ defined on the vertices of the graph) as $\hat{{\bm{f}}}={\bm{U}}^{T}{\bm{f}}$, and the inverse graph Fourier transform as ${\bm{f}}={\bm{U}}\hat{{\bm{f}}}$. Analogously to the classical Fourier transform, GFT provides a way to represent a graph signal in two domains: the vertex domain and the graph spectral domain; to filter graph signal according to smoothness; and to define the graph convolution operator, denoted as $*_{\mathcal{G}}$:

where ${\bm{g}}$ denotes the convolution kernel, and $\odot$ is the element-wise Hadamard product. If we replace the vector ${\bm{U}}^{T}{\bm{g}}$ by a diagonal matrix $\tilde{{\bm{g}}}$, then we can rewrite the Hadamard product in Eq. (8) to matrix multiplication as ${\bm{U}}\tilde{{\bm{g}}}{\bm{U}}^{T}{\bm{f}}$ (that is understood as filtering the signal ${\bm{f}}$ by the filter $\tilde{{\bm{g}}}$). Based on GFT, [1] and [38] construct convolutional neural networks (CNNs) learning on spectral domain for discrete structures such as graphs. However, there are two fundamental limitations of GFT:

• •

  High computational cost: eigendecomposition of the graph Laplacian has complexity $O(n^{3})$, and “Fourier transform” itself involves multiplying the signal with a dense matrix of eigenvectors.

• •

  The graph convolution is not localized in the vertex domain, even if the graph itself has well defined local communities.

To address these limitations, we propose a modified spectral graph network based on the MMF wavelet basis rather than the eigenbasis of the Laplacian. This has the following advantages: (i) the wavelets are generally localized in both vertex domain and frequency, (ii) the individual basis transforms are sparse, and (iii) MMF provides a computationally efficient way of decomposing graph signals into components at different granularity levels and an excellent basis for sparse approximations.

In this section, we define a convolution layer based on the wavelet bases from the MMF. This construction is inspired mainly from the GFT defined in [1] and the connection between MMF and multiresolution analysis suggested in [2].

[2] demonstrated that MMF aligns with the classical theory of multiresolution analysis (MRA), transitioning from the real line [31] to discrete spaces. In MRA of a symmetric matrix $\ignorespaces A\in\mathbb{R}^{n\times n}$ , the goal is to identify a sequence of subspaces:

The process is akin to an iterative refinement, where each subspace $\ignorespaces{\mathbb{V}}_{\ell}$ is decomposed into an orthogonal sum: $\ignorespaces{\mathbb{V}}_{\ell}={\mathbb{V}}_{\ell+1}\oplus{\mathbb{W}}_{\ell+1}$ , comprising a smoother part $\ignorespaces{\mathbb{V}}_{\ell+1}$ (the approximation space) and a rougher part $\ignorespaces{\mathbb{W}}_{\ell+1}$ (the detail space) (refer to Fig. 6). Within each subspace $\ignorespaces{\mathbb{V}}_{\ell}$ , there exists an orthonormal basis denoted by $\ignorespaces\Phi_{\ell}\triangleq\{\phi_{m}^{\ell}\}_{m}$ , where each basis function is referred to as a father wavelet. Similarly, the complementary space $\ignorespaces{\mathbb{W}}_{\ell}$ possesses an orthonormal basis denoted by $\ignorespaces\Psi_{\ell}\triangleq\{\psi_{m}^{\ell}\}_{m}$ , with each basis function termed a mother wavelet (see 11.2 for a deeper exploration of MMF’s interpretation within multiresolution analysis). Based on these wavelet bases, we can define a wavelet transform for a symmetric matrix.

In the case ${\bm{A}}$ is the normalized graph Laplacian of a graph $\mathcal{G}=(V,E)$, the wavelet transform (up to level $L$) expresses a graph signal (function over the vertex domain) $f:V\rightarrow\mathbb{R}$, without loss of generality $f\in{\mathbb{V}}_{0}$, as:

where $\alpha_{m}^{\ell}=\langle f,\psi_{m}^{\ell}\rangle$ and $\beta_{m}=\langle f,\phi_{m}^{L}\rangle$ are the wavelet coefficients. At each level, a set of coordinates ${\mathbb{T}}_{\ell}\subset{\mathbb{S}}_{\ell-1}$ are selected to be the wavelet indices, and then to be eliminated from the active set by setting ${\mathbb{S}}_{\ell}={\mathbb{S}}_{\ell-1}\setminus{\mathbb{T}}_{\ell}$ (see Section 5.1). Practically, we make the assumption that we only select $1$ wavelet index for each level that results in a single mother wavelet $\psi^{\ell}=[{\bm{A}}_{\ell}]_{i^{*},:}$ where $i^{*}$ is the selected index (see Section 5.1). We get exactly $L$ mother wavelets $\overline{\psi}=\{\psi^{1},\psi^{2},\dots,\psi^{L}\}$. On the another hand, the active rows of ${\bm{H}}={\bm{A}}_{L}$ make exactly $N-L$ father wavelets $\overline{\phi}=\{\phi^{L}_{m}={\bm{H}}_{m,:}\}_{m\in{\mathbb{S}}_{L}}$. In total, a graph of $N$ vertices has exactly $N$ wavelets (both mothers and fathers).

Analogous to the convolution based on GFT [1], each convolution layer $k=1,..,K$ of our wavelet network transforms an input vector ${\bm{f}}^{(k-1)}$ of size $\lvert V\rvert\times F_{k-1}$ into an output ${\bm{f}}^{(k)}$ of size $\lvert V\rvert\times F_{k}$ as

where ${\bm{W}}$ is our wavelet basis matrix as we concatenate $\overline{\phi}$ and $\overline{\psi}$ column-by-column, ${\bm{g}}^{(k)}_{i,j}$ is a parameter/filter in the form of a diagonal matrix learned in spectral domain similar to the filter used in the original GFT construction [1], and $\sigma$ is an element-wise linearity (e.g., ReLU, sigmoid, etc.). Each layer transforms the input features $\ignorespaces f^{(k-1)}$ into a different domain, performs filtering operations defined by the parameter $\ignorespaces g^{(k)}$ before reverting features back to the original domain. The training process is responsible for tuning the filter $\ignorespaces g^{(k)}$ to extract relevant information from the input graph signal.

We trained and evaluated our wavelet networks (WNNs) on standard graph classification benchmarks including four bioinformatics datasets: (1) MUTAG, which is a dataset of 188 mutagenic aromatic and heteroaromatic nitro compounds with 7 discrete labels [47]; (2) PTC, which consists of 344 chemical compounds with 19 discrete labels that have been tested for positive or negative toxicity in lab rats [48]; (3) PROTEINS, which contains 1,113 molecular graphs with binary labels, where nodes are secondary structure elements (SSEs) and there is an edge between two nodes if they are neighbors in the amino-acid sequence or in 3D space [49]; (4) NCI1, which has 4,110 compounds with binary labels, each screened for activity against small cell lung cancer and ovarian cancer lines [50]. Each molecule is represented by an adjacency matrix, and we represent each atomic type as a one-hot vector and use them as the node features.

We factorize all normalized graph Laplacian matrices in these datasets by MMF with $K=2$ to obtain the wavelet bases. Again, MMF wavelets are sparse and suitable for fast transform via sparse matrix multiplication. The sparsity of wavelet bases, as shown in Table 3, highlights a significant compression compared to the Fourier bases derived from the eigendecomposition of the graph Laplacian.

Our WNNs contain 6 layers of spectral convolution, 32 hidden units for each node, and are trained with 256 epochs by Adam optimization with an initial learning rate of $10^{-3}$. We follow the evaluation protocol of 10-fold cross-validation from [39]. We compare our results to several deep learning methods and popular graph kernel methods. Baseline results are taken from [46].

For graph kernel methods, we compare our model with four popular approaches: the graphlet kernel (GK), the random walk kernel (RW), the propagation kernel (PK), and the Weisfeiler-Lehman subtree kernel (WL). Each of these methods employs unique strategies for capturing graph structure and similarity. The graphlet kernel focuses on counting occurrences of small subgraphs, known as graphlets, to measure graph similarity. In contrast, the random walk kernel simulates random walks on graphs and compares the distributions of these walks to compute similarity. The propagation kernel, on the other hand, considers the propagation of labels or information through the graph to determine similarity. Lastly, the Weisfeiler-Lehman subtree kernel compares graphs based on the structural information captured by subtrees rooted at each node, iteratively refining node representations. Our results demonstrate that our method outperforms all other kernel methods on three datasets: MUTAG, PTC, and PROTEINS. However, the WL kernel achieves the best performance on the NCI1 dataset.

For deep learning methods, we compare our model with five established approaches from the literature: DGCNN, PSCN, DCNN, CCN, and IEGN. Thse deep learning methods are more closely related to our WWN, as they all leverage neural network architectures for graph analysis. Our method ranks 3rd on the MUTAG dataset, 2nd on the PTC dataset, 1st on the PROTEINS dataset, and performs the worst on the NCI1 dataset.

Our WNNs outperform 6/8, 7/8, 8/8, and 2/8 baseline methods on MUTAG, PTC, PROTEINS, and NCI1, respectively (see Table 2).

Despite these promising results, our WNN model has some limitations that we need to address in future work. One significant limitation is its performance on the NCI1 dataset, where it doese not perform as well as other methods. This suggests that our model might struggle with certain types of graphs or larger datasets. Especially considering that among the four datasets, NCI1 has the largest number of unique atom types (37, compared to less than 20 in the other three datasets) and it is also the largest dataset in terms of size. Nonetheless, further experimentation is needed to confirm the types of graph for which the model’s accuracy is suboptimal.

To further evaluate the wavelet bases returned by our learnable MMF algorithm, we construct our wavelet networks (WNNs) as in Sec. 6 and apply it to the task of node classification on two citation graphs, Cora $(N=2,708)$ and Citeseer $(N=3,312)$ [51] in which nodes and edges represent documents and citation links.

In node classification tasks, assume the number of classes is $C$, the set of labeled nodes is $V_{\text{label}}$, and we are given a normalized graph Laplacian $\tilde{L}$ and an input node feature matrix ${\bm{f}}^{(0)}$. First of all, we apply our MMF learning algorithm to factorize $\tilde{L}$ and produce our wavelet basis matrix ${\bm{W}}$. Then, we construct our wavelet network as a multi-layer CNNs with each convolution is defined as in Eq. (10) that transforms ${\bm{f}}^{(0)}$ into ${\bm{f}}^{(K)}$ after $K$ layers. The top convolution layer $K$-th returns exactly $F_{K}=C$ features and uses softmax instead of the nonlinearity $\sigma$ for each node. The loss is the cross-entropy error over all labeled nodes as:

where ${\bm{y}}_{v,c}$ is a binary indicator that is equal to $1$ if node $v$ is labeled with class $c$, and $0$ otherwise. The set of weights $\{{\bm{g}}^{(k)}\}_{k=1}^{K}$ are trained using gradient descent optimizing the loss in Eq. (11).

Each document in Cora and Citeseer has an associated feature vector (of length $1,433$ resp. $3,703$) computed from word frequencies, and is classified into one of $7$ and $6$ classes, respectively. We factorize the normalized graph Laplacian by learnable MMF with $K=16$ to obtain the wavelet bases. The resulting MMF wavelets are sparse, which makes it possible to run a fast transform on the node features by sparse matrix multiplication: only $4.69\%$ and $15.25\%$ of elements are non-zero in Citeseer and Cora, respectively. In constrast, Fourier bases given by eigendecomposition of the graph Laplacian are completely dense ($100\%$ of elements are non-zero). We evaluate our WNNs with 3 different random splits of train/validation/test: (1) $20\%$/$20\%$/$60\%$ denoted as MMF₁, (2) $40\%$/$20\%$/$40\%$ denoted as MMF₂, and (3) $60\%$/$20\%$/$20\%$ denoted as MMF₃. The WNN learns to encode the whole graph with $6$ layers of spectral convolution and $100$ hidden dimensions for each node. During training, the network is only trained to predict the node labels in the training set. Hyperparameter searching is done on the validation set. The number of epochs is $256$ and we use the Adam optimization method [52] with learning rate $\eta=10^{-3}$. We report the final test accuracy for each split in Table 4.

We compare with several traditional methods and deep learning methods including other spectral graph convolution networks such as Spectral CNN, and graph wavelet neural networks (GWNN). Baseline results are taken from [30]. Our wavelet networks perform competitively against state-of-the-art methods in the field.

We evaluate the performance of our MMF learning algorithm in comparison with the original greedy algorithm [2] and the Nyström method [61] in the task of matrix factorization on 3 datasets: (i) normalized graph Laplacian of the Karate club network ($N=34$, $E=78$) [62]; (ii) a Kronecker product matrix ($N=512$), $\mathcal{K}_{1}^{n}$, of order $n=9$, where $\mathcal{K}_{1}=((0,1),(1,1))$ is a $2\times 2$ seed matrix [63]; and (iii) normalized graph Laplacian of a Cayley tree or Bethe lattice with coordination number $z=4$ and $4$ levels of depth ($N=161$). The rotation matrix size $K$ are $8$, $16$ and $8$ for Karate, Kronecker and Cayley, respectively. Meanwhile, the original greedy MMF is limited to $K=2$ and implements an exhaustive search to find an optimal pair of indices for each rotation. For both versions of MMF, we drop $c=1$ columns after each rotation, which results in a final core size of $d_{L}=N-c\times L$. The exception is for the Kronecker matrix ($N=512$), our learning algorithm drops up to $8$ columns (for example, $L=62$ and $c=8$ results into $d_{L}=16$) to make sure that the number of learnable parameters $L\times K^{2}$ is much smaller the matrix size $N^{2}$. Our learning algorithm compresses the Kronecker matrix down to $6-7\%$ of its original size. The details of efficient training reinforcement learning with the policy networks implemented by GNNs are included in the Appendix.

For the baseline of Nyström method, we randomly select, by uniform sampling without replacement, the same number $d_{L}$ columns ${\bm{C}}$ from ${\bm{A}}$ and take out ${\bm{W}}$ as the corresponding $d_{L}\times d_{L}$ submatrix of ${\bm{A}}$. The Nyström method approximates ${\bm{A}}\approx{\bm{C}}{\bm{W}}^{\dagger}{\bm{C}}^{T}$. We measure the approximation error in Frobenius norm. Figure 5 shows our MMF learning algorithm consistently outperforms the original greedy algorithm and the Nyström baseline given the same number of active columns, $d_{L}$.