Compute Eig#

eig #

FUNCTION DESCRIPTION
eigvals

Compute eigenvalues of a square matrix.
eigvecs

Compute eigenvalues and eigenvectors of a square matrix.
gen_eig_chol

Solve a symmetric generalized eigenvalue problem using Cholesky

Functions#

eigvals(X: torch.Tensor, n_eigen: int = 0, symmetric: bool = True, sort_descending: bool = True, UPLO: Literal['U', 'L'] = 'U') -> torch.Tensor #

Compute eigenvalues of a square matrix.

PARAMETER DESCRIPTION
X

Square matrix of shape (n_samples, n_samples). For symmetric equal to True, matrix should be Hermitian.

TYPE: Tensor

n_eigen

Number of eigenvalues to return. If 0, returns all eigenvalues (n_samples, ).

TYPE: int, optional, by default 0 DEFAULT: 0

symmetric

Eigendecomposition mode. Set to True for symmetric or Hermitian matrices.

TYPE: bool, optional, by default True DEFAULT: True

sort_descending

Whether to sort eigenvalues in descending order.

TYPE: bool, optional, by default True DEFAULT: True

UPLO

Whether to use upper ("U") or lower ("L") triangular part for symmetric matrices.

TYPE: Literal["U", "L"], optional, by default "U" DEFAULT: 'U'

RETURNS DESCRIPTION
Tensor

Eigenvalues of shape (n_eigen,) or (n_samples, ) if number of eigenvalues (n_eigen) is set to 0.
RAISES DESCRIPTION
ValueError

If n_eigen is negative or exceeds matrix size.
RuntimeError

If matrix is not square or not 2D.

Examples:

Compute all eigenvalues of a symmetric matrix:

>>> import torch
>>> X = torch.tensor([[4.0, 2.0], [2.0, 3.0]])
>>> eigenvalues = eigvals(X)

Get top 2 eigenvalues:

>>> top_eigenvalues = eigvals(X, n_eigen=2)

Eigenvalues of a nonsymmetric matrix:

>>> X_nonsym = torch.tensor([[1.0, 2.0], [3.0, 4.0]])
>>> eigenvalues = eigvals(X_nonsym, symmetric=False)
Source code in spectre/compute/eig.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def eigvals(
    X: torch.Tensor,
    n_eigen: int = 0,
    symmetric: bool = True,
    sort_descending: bool = True,
    UPLO: Literal["U", "L"] = "U",
) -> torch.Tensor:
    """
    Compute eigenvalues of a square matrix.

    Parameters
    ----------
    X : torch.Tensor
        Square matrix of shape (n_samples, n_samples). For symmetric equal to
        True, matrix should be Hermitian.

    n_eigen : int, optional, by default 0
        Number of eigenvalues to return. If 0, returns all eigenvalues
        (n_samples, ).

    symmetric : bool, optional, by default True
        Eigendecomposition mode. Set to True for symmetric or Hermitian
        matrices.

    sort_descending : bool, optional, by default True
        Whether to sort eigenvalues in descending order.

    UPLO : Literal["U", "L"], optional, by default "U"
        Whether to use upper ("U") or lower ("L") triangular part for symmetric
        matrices.

    Returns
    -------
    torch.Tensor
        Eigenvalues of shape (n_eigen,) or (n_samples, ) if number of
        eigenvalues (n_eigen) is set to 0.

    Raises
    ------
    ValueError
        If n_eigen is negative or exceeds matrix size.

    RuntimeError
        If matrix is not square or not 2D.

    Examples
    --------
    Compute all eigenvalues of a symmetric matrix:

    >>> import torch
    >>> X = torch.tensor([[4.0, 2.0], [2.0, 3.0]])
    >>> eigenvalues = eigvals(X)

    Get top 2 eigenvalues:

    >>> top_eigenvalues = eigvals(X, n_eigen=2)

    Eigenvalues of a nonsymmetric matrix:

    >>> X_nonsym = torch.tensor([[1.0, 2.0], [3.0, 4.0]])
    >>> eigenvalues = eigvals(X_nonsym, symmetric=False)
    """
    check_2d(X)
    check_square_shape(X)

    n_samples = X.shape[0]

    if not isinstance(n_eigen, int):
        raise TypeError("Number of eigenvalues must be an integer.")
    check_in_interval(n_eigen, interval=f"[0,{n_samples}]")

    if symmetric:
        eigval = torch.linalg.eigvalsh(X, UPLO=UPLO)
    else:
        eigval = torch.linalg.eigvals(X)
        eigval = torch.real(eigval)

    eigval = eigval[torch.argsort(eigval, descending=sort_descending)]

    if n_eigen == 0:
        return eigval

    return eigval[:n_eigen]

eigvecs(X: torch.Tensor, n_eigen: int = 0, symmetric: bool = True, sort_descending: bool = True, UPLO: Literal['U', 'L'] = 'U') -> tuple[torch.Tensor, torch.Tensor] #

Compute eigenvalues and eigenvectors of a square matrix.

PARAMETER DESCRIPTION
X

Square matrix of shape (n_samples, n_samples). For symmetric equal to True, matrix should be Hermitian.

TYPE: Tensor

n_eigen

Number of eigenvalues to return. If 0, returns all eigenvalues (n_samples, ).

TYPE: int, optional, by default 0 DEFAULT: 0

symmetric

Eigendecomposition mode. Set to True for symmetric or Hermitian matrices.

TYPE: bool, optional, by default True DEFAULT: True

sort_descending

Whether to sort eigenvalues in descending order.

TYPE: bool, optional, by default True DEFAULT: True

UPLO

Whether to use upper ("U") or lower ("L") triangular part for symmetric matrices.

TYPE: Literal["U", "L"], optional, by default "U" DEFAULT: 'U'

RETURNS DESCRIPTION
tuple[Tensor, Tensor]

Eigenvalues of shape (n_eigen, ) and eigenvectors of shape (n_samples, n_eigen).
RAISES DESCRIPTION
ValueError

If n_eigen is negative or exceeds matrix size.
RuntimeError

If matrix is not square or not 2D.
Notes

Eigenvectors are normalized to canonical representation where the sum of each eigenvector is positive. This normalization ensures consistency across different runs and is particularly useful for spectral methods where eigenvector signs are arbitrary. The sign is chosen such that $\sum_i v_i

0$ for each eigenvector \(v\).

Examples:

Compute eigenvalues and eigenvectors of a symmetric matrix:

>>> import torch
>>> X = torch.tensor([[4.0, 2.0], [2.0, 3.0]])
>>> eigenvalues, eigenvectors = eigvecs(X)

Get top 2 eigenvalues and corresponding eigenvectors:

>>> eigenvals, eigenvecs_subset = eigvecs(X, n_eigen=2)

Eigendecomposition of a nonsymmetric matrix:

>>> X_nonsym = torch.tensor([[1.0, 2.0], [3.0, 4.0]])
>>> eigenvalues, eigenvectors = eigvecs(X_nonsym, symmetric=False)
Source code in spectre/compute/eig.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def eigvecs(
    X: torch.Tensor,
    n_eigen: int = 0,
    symmetric: bool = True,
    sort_descending: bool = True,
    UPLO: Literal["U", "L"] = "U",
) -> tuple[torch.Tensor, torch.Tensor]:
    """
    Compute eigenvalues and eigenvectors of a square matrix.

    Parameters
    ----------
    X : torch.Tensor
        Square matrix of shape (n_samples, n_samples). For symmetric equal to
        True, matrix should be Hermitian.

    n_eigen : int, optional, by default 0
        Number of eigenvalues to return. If 0, returns all eigenvalues
        (n_samples, ).

    symmetric : bool, optional, by default True
        Eigendecomposition mode. Set to True for symmetric or Hermitian
        matrices.

    sort_descending : bool, optional, by default True
        Whether to sort eigenvalues in descending order.

    UPLO : Literal["U", "L"], optional, by default "U"
        Whether to use upper ("U") or lower ("L") triangular part for symmetric
        matrices.

    Returns
    -------
    tuple[torch.Tensor, torch.Tensor]
        Eigenvalues of shape (n_eigen, ) and eigenvectors of shape (n_samples,
        n_eigen).

    Raises
    ------
    ValueError
        If n_eigen is negative or exceeds matrix size.

    RuntimeError
        If matrix is not square or not 2D.

    Notes
    -----
    Eigenvectors are normalized to canonical representation where the sum of
    each eigenvector is positive. This normalization ensures consistency across
    different runs and is particularly useful for spectral methods where
    eigenvector signs are arbitrary. The sign is chosen such that $\\sum_i v_i
    > 0$ for each eigenvector $v$.

    Examples
    --------
    Compute eigenvalues and eigenvectors of a symmetric matrix:

    >>> import torch
    >>> X = torch.tensor([[4.0, 2.0], [2.0, 3.0]])
    >>> eigenvalues, eigenvectors = eigvecs(X)

    Get top 2 eigenvalues and corresponding eigenvectors:

    >>> eigenvals, eigenvecs_subset = eigvecs(X, n_eigen=2)

    Eigendecomposition of a nonsymmetric matrix:

    >>> X_nonsym = torch.tensor([[1.0, 2.0], [3.0, 4.0]])
    >>> eigenvalues, eigenvectors = eigvecs(X_nonsym, symmetric=False)
    """
    check_2d(X)
    check_square_shape(X)

    n_samples = X.shape[0]

    if not isinstance(n_eigen, int):
        raise TypeError("Number of eigenvalues must be an integer.")
    check_in_interval(n_eigen, interval=f"[0,{n_samples}]")

    if symmetric:
        eigval, eigvec = torch.linalg.eigh(X, UPLO=UPLO)
    else:
        eigval, eigvec = torch.linalg.eig(X)
        eigval = torch.real(eigval)
        eigvec = torch.real(eigvec)

    idx = torch.argsort(eigval, descending=sort_descending)
    eigval = eigval[idx]
    eigvec = eigvec[:, idx]

    # Convert eigenvectors to canonical representation (sum > 0).
    eigvec[:, eigvec.sum(dim=0) < 0] *= -1

    if n_eigen == 0:
        return eigval, eigvec

    return eigval[:n_eigen], eigvec[:, :n_eigen]

gen_eig_chol(C: torch.Tensor, Q: torch.Tensor, n_eigen: int = 0, sort_descending: bool = True, eps: float = torch.finfo(torch.float32).eps) -> tuple[torch.Tensor, torch.Tensor] #

Solve a symmetric generalized eigenvalue problem using Cholesky decomposition [1].

Converts the generalized eigenvalue problem \(C v = w Q v\) to a regular symmetric eigenvalue problem using Cholesky factorization of \(Q\).

Algorithm:

  • Perform Cholesky decomposition \(Q = L R\), where \(R=L^\top\).
  • Calculate \(C' = L^{-1} C R^{-1}\).
  • Solve a standard eigenvalue problem \(C' v = w v'\).
  • Recover original eigenvectors and eigenvalues \(v = R^{-1} v'\).
PARAMETER DESCRIPTION
C

Symmetric matrix of shape (n_samples, n_samples).

TYPE: Tensor

Q

Symmetric positive semi-definite matrix of shape (n_samples, n_samples).

TYPE: Tensor

n_eigen

Number of eigenvalues to return. If 0, returns all eigenvalues (n_samples, ).

TYPE: int, optional, by default 0 DEFAULT: 0

sort_descending

Whether to sort eigenvalues in descending order.

TYPE: bool, optional, by default True DEFAULT: True

eps

Small value added to diagonal of Q for numerical stability.

TYPE: float, optional, by default torch.finfo(torch.float32).eps DEFAULT: eps

RETURNS DESCRIPTION
tuple[Tensor, Tensor]

Eigenvalues in descending order and corresponding eigenvectors with shape (n_samples, ) and (n_samples, n_samples), respectively.
RAISES DESCRIPTION
ValueError

If Q matrix is not positive semi-definite.
Source code in spectre/compute/eig.py 
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
def gen_eig_chol(
    C: torch.Tensor,
    Q: torch.Tensor,
    n_eigen: int = 0,
    sort_descending: bool = True,
    eps: float = torch.finfo(torch.float32).eps,
) -> tuple[torch.Tensor, torch.Tensor]:
    """
    Solve a symmetric generalized eigenvalue problem using Cholesky
    decomposition [@chen2019nonlinear].

    Converts the generalized eigenvalue problem $C v = w Q v$ to a regular
    symmetric eigenvalue problem using Cholesky factorization of $Q$.

    Algorithm:

    - Perform Cholesky decomposition $Q = L R$, where $R=L^\\top$.
    - Calculate $C' = L^{-1} C R^{-1}$.
    - Solve a standard eigenvalue problem $C' v = w v'$.
    - Recover original eigenvectors and eigenvalues $v = R^{-1} v'$.

    Parameters
    ----------
    C : torch.Tensor
        Symmetric matrix of shape (n_samples, n_samples).

    Q : torch.Tensor
        Symmetric positive semi-definite matrix of shape (n_samples,
        n_samples).

    n_eigen : int, optional, by default 0
        Number of eigenvalues to return. If 0, returns all eigenvalues
        (n_samples, ).

    sort_descending : bool, optional, by default True
        Whether to sort eigenvalues in descending order.

    eps : float, optional, by default torch.finfo(torch.float32).eps
        Small value added to diagonal of Q for numerical stability.

    Returns
    -------
    tuple[torch.Tensor, torch.Tensor]
        Eigenvalues in descending order and corresponding eigenvectors with
        shape (n_samples, ) and (n_samples, n_samples), respectively.

    Raises
    ------
    ValueError
        If Q matrix is not positive semi-definite.
    """
    check_2d(C)
    check_square_shape(C)
    check_2d(Q)
    check_square_shape(Q)

    n_samples = C.shape[0]

    if not isinstance(n_eigen, int):
        raise TypeError("Number of eigenvalues must be an integer.")
    check_in_interval(n_eigen, interval=f"[0,{n_samples}]")

    # `torch.linalg.cholesky` checks for Hermitian matrix automatically and
    # throws runtime error if violated.
    try:
        Q = Q + eps * torch.eye(Q.shape[0], device=Q.device)
        L = torch.linalg.cholesky(Q)
        L_inv = torch.linalg.inv(L)
        R_inv = torch.linalg.inv(L.t())

        C_tilde = torch.matmul(L_inv, torch.matmul(C, R_inv))

        # Regular symmetric eigenvalue problem.
        w, v_tilde = torch.linalg.eigh(C_tilde)

        # Recover original eigenvalues and eigenvectors.
        v = torch.matmul(R_inv, v_tilde)
    except RuntimeError as e:
        raise ValueError(
            f"Q matrix is not positive semi-definite. Cholesky decomposition failed: {e}"
        ) from e

    # Reordering to non-ascending.
    if sort_descending:
        w = torch.flip(w, dims=[0])
        v = torch.flip(v, dims=[1])

    if n_eigen == 0:
        return w, v

    return w[:n_eigen], v[:, :n_eigen]