Coverage for src/cvx/linalg/decomposition/power_iteration.py: 100%

1"""Dominant eigenpair estimation via power iteration.""" 

2 

3from __future__ import annotations 

4 

5from collections.abc import Callable 

6 

7import numpy as np 

8 

9from ..core.exceptions import NonSquareMatrixError, NotAMatrixError 

10from ..core.types import Matrix, SupportsMatvec, Vector 

11 

12_CALLABLE_NEEDS_N_MESSAGE = "power_iteration needs `n` when `operator` is a bare callable" 

13 

14 

15def power_iteration( 

16 operator: Matrix | SupportsMatvec | Callable[[Vector], Vector], 

17 *, 

18 n: int | None = None, 

19 n_iter: int = 1000, 

20 tol: float = 1e-9, 

21 seed: int | None = None, 

22) -> tuple[float, Vector]: 

23 """Estimate the dominant eigenpair of a real symmetric operator. 

24 

25 Repeatedly applies *operator* to a random unit vector, renormalizing each 

26 step, until the Rayleigh-quotient eigenvalue estimate stops changing. The 

27 iterate converges to the eigenvector whose eigenvalue is largest in 

28 magnitude; the matching eigenvalue (returned with its sign) comes from the 

29 Rayleigh quotient ``v.T @ (operator @ v)``. 

30 

31 *operator* may be given three ways, so the leading eigenvalue can be 

32 estimated **matrix-free** (e.g. the Lipschitz constant of a gradient step): 

33 

34 * a dense symmetric ``(n, n)`` array; 

35 * any :class:`~cvx.linalg.SymmetricOperator` (its :meth:`matvec` and ``n`` 

36 drive the iteration, so no ``n x n`` matrix is formed); or 

37 * a callable ``v -> A @ v`` together with the dimension *n*. 

38 

39 Each iteration costs a single application, so this is far cheaper than a full 

40 :func:`~cvx.linalg.eigh` when only the leading eigenpair is needed and there 

41 is a clear spectral gap. Symmetry is assumed; only the result's 

42 interpretation as an eigenpair relies on it. Like :func:`~cvx.linalg.svd`, 

43 this is a raw primitive and is **not** NaN-aware. 

44 

45 Args: 

46 operator: A symmetric matrix, a :class:`~cvx.linalg.SymmetricOperator`, 

47 or a callable applying it to a vector. 

48 n: Dimension of the operator. Required when *operator* is a bare callable; 

49 ignored otherwise (taken from the array shape or the operator's ``n``). 

50 n_iter: Maximum number of iterations. Defaults to 1000. 

51 tol: Convergence tolerance on the relative change of the eigenvalue 

52 estimate between iterations. Defaults to ``1e-9``. 

53 seed: Seed for the random starting vector, for reproducibility. If 

54 ``None``, uses the current NumPy random state. 

55 

56 Returns: 

57 Tuple ``(eigenvalue, eigenvector)`` where *eigenvalue* is the signed 

58 dominant eigenvalue estimate (a ``float``) and *eigenvector* is the 

59 corresponding unit-norm eigenvector. The eigenvector sign is arbitrary. 

60 

61 Raises: 

62 NotAMatrixError: If *operator* is an array that is not 2-D. 

63 NonSquareMatrixError: If *operator* is an array that is not square. 

64 ValueError: If *operator* is a bare callable and *n* is not given. 

65 

66 Example: 

67 >>> import numpy as np 

68 >>> from cvx.linalg import power_iteration 

69 >>> matrix = np.diag([3.0, 2.0, 1.0]) 

70 >>> eigenvalue, eigenvector = power_iteration(matrix, seed=0) 

71 >>> bool(np.isclose(eigenvalue, 3.0)) 

72 True 

73 >>> bool(np.isclose(abs(eigenvector[0]), 1.0)) 

74 True 

75 

76 It runs matrix-free on a :class:`~cvx.linalg.SymmetricOperator`, so the 

77 leading eigenvalue of ``M.T @ M`` needs no ``n x n`` matrix: 

78 

79 >>> from cvx.linalg import GramOperator 

80 >>> rng = np.random.default_rng(0) 

81 >>> M = rng.standard_normal((20, 5)) 

82 >>> lam, _ = power_iteration(GramOperator(M), seed=0) 

83 >>> bool(np.isclose(lam, np.linalg.eigvalsh(M.T @ M)[-1])) 

84 True 

85 """ 

86 apply, dim = _resolve(operator, n) 

87 

88 rng = np.random.default_rng(seed) 

89 v = rng.standard_normal((dim,)) 

90 v = v / np.linalg.norm(v) 

91 

92 eigenvalue = float(v @ apply(v)) 

93 for _ in range(n_iter): 

94 w = apply(v) 

95 norm = float(np.linalg.norm(w)) 

96 if norm < 1e-15: 

97 # operator annihilates the iterate: the dominant eigenvalue is zero. 

98 return 0.0, v 

99 v = w / norm 

100 new_eigenvalue = float(v @ apply(v)) 

101 if abs(new_eigenvalue - eigenvalue) <= tol * max(1.0, abs(new_eigenvalue)): 

102 return new_eigenvalue, v 

103 eigenvalue = new_eigenvalue 

104 

105 return eigenvalue, v 

106 

107 

108def _resolve( 

109 operator: Matrix | SupportsMatvec | Callable[[Vector], Vector], 

110 n: int | None, 

111) -> tuple[Callable[[Vector], Vector], int]: 

112 """Return ``(apply, dimension)`` for an array, a SymmetricOperator, or a callable. 

113 

114 The operator is matched structurally against :class:`~cvx.linalg.core.types.SupportsMatvec` 

115 (a lower layer than :mod:`~cvx.linalg.operators`), so this needs no import of the 

116 operator backends and closes no import cycle. 

117 """ 

118 if isinstance(operator, SupportsMatvec): 

119 return operator.matvec, operator.n 

120 if not isinstance(operator, np.ndarray) and callable(operator): 

121 if n is None: 

122 raise ValueError(_CALLABLE_NEEDS_N_MESSAGE) 

123 return operator, int(n) 

124 return _array_apply(operator) 

125 

126 

127def _array_apply( 

128 operator: Matrix | SupportsMatvec | Callable[[Vector], Vector], 

129) -> tuple[Callable[[Vector], Vector], int]: 

130 """Return ``(v -> A @ v, n)`` for a dense square array, validating its shape.""" 

131 array = np.asarray(operator) 

132 if array.ndim != 2: 

133 raise NotAMatrixError(array.ndim, func="power_iteration") 

134 rows, cols = array.shape 

135 if rows != cols: 

136 raise NonSquareMatrixError(rows, cols) 

137 return lambda v: array @ v, cols