|
| 1 | +""" |
| 2 | +Lanczos Method for Finding Eigenvalues and Eigenvectors of a Graph. |
| 3 | +
|
| 4 | +This module demonstrates the Lanczos method to approximate the largest eigenvalues |
| 5 | +and corresponding eigenvectors of a symmetric matrix represented as a graph's |
| 6 | +adjacency list. The method efficiently handles large, sparse matrices by converting |
| 7 | +the graph to a tridiagonal matrix, whose eigenvalues and eigenvectors are then |
| 8 | +computed. |
| 9 | +
|
| 10 | +Key Functions: |
| 11 | +- `find_lanczos_eigenvectors`: Computes the k largest eigenvalues and vectors. |
| 12 | +- `lanczos_iteration`: Constructs the tridiagonal matrix and orthonormal basis vectors. |
| 13 | +- `multiply_matrix_vector`: Multiplies an adjacency list graph with a vector. |
| 14 | +
|
| 15 | +Complexity: |
| 16 | +- Time: O(k * n), where k is the number of eigenvalues and n is the matrix size. |
| 17 | +- Space: O(n), due to sparse representation and tridiagonal matrix structure. |
| 18 | +
|
| 19 | +Further Reading: |
| 20 | +- Lanczos Algorithm: https://en.wikipedia.org/wiki/Lanczos_algorithm |
| 21 | +- Eigenvector Centrality: https://en.wikipedia.org/wiki/Eigenvector_centrality |
| 22 | +
|
| 23 | +Example Usage: |
| 24 | +Given a graph represented by an adjacency list, the `find_lanczos_eigenvectors` |
| 25 | +function returns the largest eigenvalues and eigenvectors. This can be used to |
| 26 | +analyze graph centrality. |
| 27 | +""" |
| 28 | + |
| 29 | +import numpy as np |
| 30 | + |
| 31 | + |
| 32 | +def validate_adjacency_list(graph: list[list[int | None]]) -> None: |
| 33 | + """Validates the adjacency list format for the graph. |
| 34 | +
|
| 35 | + Args: |
| 36 | + graph: A list of lists where each sublist contains the neighbors of a node. |
| 37 | +
|
| 38 | + Raises: |
| 39 | + ValueError: If the graph is not a list of lists, or if any node has |
| 40 | + invalid neighbors (e.g., out-of-range or non-integer values). |
| 41 | +
|
| 42 | + >>> validate_adjacency_list([[1, 2], [0], [0, 1]]) |
| 43 | + >>> validate_adjacency_list([[]]) # No neighbors, valid case |
| 44 | + >>> validate_adjacency_list([[1], [2], [-1]]) # Invalid neighbor |
| 45 | + Traceback (most recent call last): |
| 46 | + ... |
| 47 | + ValueError: Invalid neighbor -1 in node 2 adjacency list. |
| 48 | + """ |
| 49 | + if not isinstance(graph, list): |
| 50 | + raise ValueError("Graph should be a list of lists.") |
| 51 | + |
| 52 | + for node_index, neighbors in enumerate(graph): |
| 53 | + if not isinstance(neighbors, list): |
| 54 | + no_neighbors_message: str = ( |
| 55 | + f"Node {node_index} should have a list of neighbors." |
| 56 | + ) |
| 57 | + raise ValueError(no_neighbors_message) |
| 58 | + for neighbor_index in neighbors: |
| 59 | + if ( |
| 60 | + not isinstance(neighbor_index, int) |
| 61 | + or neighbor_index < 0 |
| 62 | + or neighbor_index >= len(graph) |
| 63 | + ): |
| 64 | + invalid_neighbor_message: str = ( |
| 65 | + f"Invalid neighbor {neighbor_index} in node {node_index} " |
| 66 | + f"adjacency list." |
| 67 | + ) |
| 68 | + raise ValueError(invalid_neighbor_message) |
| 69 | + |
| 70 | + |
| 71 | +def lanczos_iteration( |
| 72 | + graph: list[list[int | None]], num_eigenvectors: int |
| 73 | +) -> tuple[np.ndarray, np.ndarray]: |
| 74 | + """Constructs the tridiagonal matrix and orthonormal basis vectors using the |
| 75 | + Lanczos method. |
| 76 | +
|
| 77 | + Args: |
| 78 | + graph: The graph represented as a list of adjacency lists. |
| 79 | + num_eigenvectors: The number of largest eigenvalues and eigenvectors |
| 80 | + to approximate. |
| 81 | +
|
| 82 | + Returns: |
| 83 | + A tuple containing: |
| 84 | + - tridiagonal_matrix: A (num_eigenvectors x num_eigenvectors) symmetric |
| 85 | + matrix. |
| 86 | + - orthonormal_basis: A (num_nodes x num_eigenvectors) matrix of orthonormal |
| 87 | + basis vectors. |
| 88 | +
|
| 89 | + Raises: |
| 90 | + ValueError: If num_eigenvectors is less than 1 or greater than the number of |
| 91 | + nodes. |
| 92 | +
|
| 93 | + >>> graph = [[1, 2], [0, 2], [0, 1]] |
| 94 | + >>> T, Q = lanczos_iteration(graph, 2) |
| 95 | + >>> T.shape == (2, 2) and Q.shape == (3, 2) |
| 96 | + True |
| 97 | + """ |
| 98 | + num_nodes: int = len(graph) |
| 99 | + if not (1 <= num_eigenvectors <= num_nodes): |
| 100 | + raise ValueError( |
| 101 | + "Number of eigenvectors must be between 1 and the number of " |
| 102 | + "nodes in the graph." |
| 103 | + ) |
| 104 | + |
| 105 | + orthonormal_basis: np.ndarray = np.zeros((num_nodes, num_eigenvectors)) |
| 106 | + tridiagonal_matrix: np.ndarray = np.zeros((num_eigenvectors, num_eigenvectors)) |
| 107 | + |
| 108 | + rng = np.random.default_rng() |
| 109 | + initial_vector: np.ndarray = rng.random(num_nodes) |
| 110 | + initial_vector /= np.sqrt(np.dot(initial_vector, initial_vector)) |
| 111 | + orthonormal_basis[:, 0] = initial_vector |
| 112 | + |
| 113 | + prev_beta: float = 0.0 |
| 114 | + for iter_index in range(num_eigenvectors): |
| 115 | + result_vector: np.ndarray = multiply_matrix_vector( |
| 116 | + graph, orthonormal_basis[:, iter_index] |
| 117 | + ) |
| 118 | + if iter_index > 0: |
| 119 | + result_vector -= prev_beta * orthonormal_basis[:, iter_index - 1] |
| 120 | + alpha_value: float = np.dot(orthonormal_basis[:, iter_index], result_vector) |
| 121 | + result_vector -= alpha_value * orthonormal_basis[:, iter_index] |
| 122 | + |
| 123 | + prev_beta = np.sqrt(np.dot(result_vector, result_vector)) |
| 124 | + if iter_index < num_eigenvectors - 1 and prev_beta > 1e-10: |
| 125 | + orthonormal_basis[:, iter_index + 1] = result_vector / prev_beta |
| 126 | + tridiagonal_matrix[iter_index, iter_index] = alpha_value |
| 127 | + if iter_index < num_eigenvectors - 1: |
| 128 | + tridiagonal_matrix[iter_index, iter_index + 1] = prev_beta |
| 129 | + tridiagonal_matrix[iter_index + 1, iter_index] = prev_beta |
| 130 | + return tridiagonal_matrix, orthonormal_basis |
| 131 | + |
| 132 | + |
| 133 | +def multiply_matrix_vector( |
| 134 | + graph: list[list[int | None]], vector: np.ndarray |
| 135 | +) -> np.ndarray: |
| 136 | + """Performs multiplication of a graph's adjacency list representation with a vector. |
| 137 | +
|
| 138 | + Args: |
| 139 | + graph: The adjacency list of the graph. |
| 140 | + vector: A 1D numpy array representing the vector to multiply. |
| 141 | +
|
| 142 | + Returns: |
| 143 | + A numpy array representing the product of the adjacency list and the vector. |
| 144 | +
|
| 145 | + Raises: |
| 146 | + ValueError: If the vector's length does not match the number of nodes in the |
| 147 | + graph. |
| 148 | +
|
| 149 | + >>> multiply_matrix_vector([[1, 2], [0, 2], [0, 1]], np.array([1, 1, 1])) |
| 150 | + array([2., 2., 2.]) |
| 151 | + >>> multiply_matrix_vector([[1, 2], [0, 2], [0, 1]], np.array([0, 1, 0])) |
| 152 | + array([1., 0., 1.]) |
| 153 | + """ |
| 154 | + num_nodes: int = len(graph) |
| 155 | + if vector.shape[0] != num_nodes: |
| 156 | + raise ValueError("Vector length must match the number of nodes in the graph.") |
| 157 | + |
| 158 | + result: np.ndarray = np.zeros(num_nodes) |
| 159 | + for node_index, neighbors in enumerate(graph): |
| 160 | + for neighbor_index in neighbors: |
| 161 | + result[node_index] += vector[neighbor_index] |
| 162 | + return result |
| 163 | + |
| 164 | + |
| 165 | +def find_lanczos_eigenvectors( |
| 166 | + graph: list[list[int | None]], num_eigenvectors: int |
| 167 | +) -> tuple[np.ndarray, np.ndarray]: |
| 168 | + """Computes the largest eigenvalues and their corresponding eigenvectors using the |
| 169 | + Lanczos method. |
| 170 | +
|
| 171 | + Args: |
| 172 | + graph: The graph as a list of adjacency lists. |
| 173 | + num_eigenvectors: Number of largest eigenvalues and eigenvectors to compute. |
| 174 | +
|
| 175 | + Returns: |
| 176 | + A tuple containing: |
| 177 | + - eigenvalues: 1D array of the largest eigenvalues in descending order. |
| 178 | + - eigenvectors: 2D array where each column is an eigenvector corresponding |
| 179 | + to an eigenvalue. |
| 180 | +
|
| 181 | + Raises: |
| 182 | + ValueError: If the graph format is invalid or num_eigenvectors is out of bounds. |
| 183 | +
|
| 184 | + >>> eigenvalues, eigenvectors = find_lanczos_eigenvectors( |
| 185 | + ... [[1, 2], [0, 2], [0, 1]], 2 |
| 186 | + ... ) |
| 187 | + >>> len(eigenvalues) == 2 and eigenvectors.shape[1] == 2 |
| 188 | + True |
| 189 | + """ |
| 190 | + validate_adjacency_list(graph) |
| 191 | + tridiagonal_matrix, orthonormal_basis = lanczos_iteration(graph, num_eigenvectors) |
| 192 | + eigenvalues, eigenvectors = np.linalg.eigh(tridiagonal_matrix) |
| 193 | + return eigenvalues[::-1], np.dot(orthonormal_basis, eigenvectors[:, ::-1]) |
| 194 | + |
| 195 | + |
| 196 | +def main() -> None: |
| 197 | + """ |
| 198 | + Main driver function for testing the implementation with doctests. |
| 199 | + """ |
| 200 | + import doctest |
| 201 | + |
| 202 | + doctest.testmod() |
| 203 | + |
| 204 | + |
| 205 | +if __name__ == "__main__": |
| 206 | + main() |
0 commit comments