33
44
55def inner_outer_eigs (eigs , rho ):
6+ """
7+ Splits the eigenvalues into inner eigenvalues and the outer eigenvalue based on the sign of rho.
8+
9+ Parameters:
10+ eigs (np.ndarray or scipy.sparse.spmatrix): Array of eigenvalues, assumed to be sorted.
11+ rho (float): Scalar parameter appearing in the secular function (please refer to the documentation for more detailed info).
12+
13+ Returns:
14+ tuple: A tuple (inner_eigs, outer_eig) where inner_eigs is an array of eigenvalues and outer_eig is a scalar.
15+ If rho > 0, the last eigenvalue is considered outer due to the interlacing property; otherwise, the first is.
16+ """
617 inner_eigs = eigs [:- 1 ] if rho > 0 else eigs [1 :]
718 outer_eig = eigs [- 1 ] if rho > 0 else eigs [0 ]
819 return inner_eigs , outer_eig
920
1021
1122def return_secular_f (rho , d , v ):
1223 """
13- This returns f as a callable object (function of lambda). f is built using rho, d, v.
14- This function passes the tests in test.py and is likely implemented correctly.
24+ Constructs the secular function for a rank-one update to a diagonal matrix.
25+
26+ Parameters:
27+ rho (float): Scalar from the rank-one matrix update.
28+ d (np.ndarray or scipy.sparse.spmatrix): 1D array or sparse vector of diagonal entries.
29+ v (np.ndarray or scipy.sparse.spmatrix): 1D array or sparse vector used in the rank-one update.
30+
31+ Returns:
32+ callable:
33+ f(lambda_: float) -> float: The secular function evaluated at lambda_.
1534 """
1635
1736 def f (lambda_ ):
@@ -25,27 +44,53 @@ def f(lambda_):
2544
2645def secular_function (mu , rho , d , v2 , i ):
2746 """
28- Needed to compute the zeros of the secular function in the i-th subinterval;
47+ Evaluates the secular function at a given point in the i-th subinterval.
48+
49+ Parameters:
50+ mu (float): Point at which to evaluate the secular function.
51+ rho (float): Scalar from the rank-one matrix update.
52+ d (np.ndarray or scipy.sparse.spmatrix): 1D array or sparse vector of diagonal entries.
53+ v2 (np.ndarray or scipy.sparse.spmatrix): Elementwise square of the update vector v (i.e., v ** 2).
54+ i (int): Index of the subinterval in which mu lies.
55+
56+ Returns:
57+ float: The value of the secular function at mu.
2958 """
3059 psi1 , _ , psi2 , _ = compute_psi_s (mu , rho , d , v2 , i )
3160 return 1 + psi1 + psi2
3261
3362
3463def check_is_root (f , x , tol = 1e-7 ):
3564 """
36- Usually the values of f(found_eigenvalue) are around 1e-10, so we cannot be *too* restrictive with the threshold.
37- For instance, using np.isclose, sometimes the checks are not passed, even though f(found_eig) is very small in
38- absolute value and we are indeed very close to one.
39- That is the reason for defining a helper function.
65+ Determines whether x is a root of the function f within a given numerical tolerance.
66+ Written because np.isclose is too restrictive even in cases in which we are indeed close to an eigenvalue.
67+
68+ Parameters:
69+ f (callable): Function to evaluate.
70+ x (float): Point to test as a root.
71+ tol (float): Absolute tolerance for considering f(x) close to zero.
72+
73+ Returns:
74+ bool: True if |f(x)| < tol, indicating x is a root within the specified tolerance.
4075 """
4176 return np .abs (f (x )) < tol
4277
4378
4479def bisection (f , a , b , tol , max_iter ):
4580 """
46- In the main method, we used a slightly tweaked form of bisection for the eig. in the outer interval (i.e. lambda_0 if rho < 0,
47- else lambda_{n-1}.
48- This helper function implements standard bisection and it is used in the function compute_outer_zero.
81+ Standard bisection method to find a root of the function f in the interval [a, b].
82+
83+ This implementation is used in `compute_outer_zero` to locate the outer eigenvalue.
84+
85+ Parameters:
86+ f (callable): A continuous function for which f(a) * f(b) < 0.
87+ a (float): Left endpoint of the interval.
88+ b (float): Right endpoint of the interval.
89+ tol (float): Tolerance for convergence. The method stops when the interval is smaller than tol or when f(c) is sufficiently small.
90+ max_iter (int): Maximum number of iterations before stopping.
91+
92+ Returns:
93+ float: Approximation of the root within the specified tolerance.
4994 """
5095 iter_count = 0
5196
@@ -65,17 +110,27 @@ def bisection(f, a, b, tol, max_iter):
65110
66111def compute_outer_zero (f , rho , interval_end , v , tol = 1e-12 , max_iter = 2000 ):
67112 """
68- Function to compute the outer eigenvalue (lambda[0] if rho < 0, else lambda[n-1]).
69- What it does is the following:
70- 1) depending on rho, understand whether we should look for it in (d[n-1],+ \infty) or (-\infty, d[0])
71- 2) use bisection as follows:
72- 2a) Fix the bisection interval. Notice that (assuming rho > 0, else it is equivalent but specular)
73- f(d[n-1]+\epsilon)\a pprox -\infty, in particular f(d[n-1]+\epsilon)<0.
74- So we just need to find r\in\mathbb{R} such that f(d[n-1]+r)>0, and we can regularly use bisection.
75- 2b) To find that value of r, we just add arbitrary values to d[n-1] until the condition is satisfied.
76- 2c) Bisection is then used
77-
78- Notice that this function passes all the tests, and I assume it is implemented correctly.
113+ Computes the outer eigenvalue (lambda[0] if rho < 0, lambda[n-1] if rho > 0) of a rank-one modified diagonal matrix.
114+
115+ The secular function behaves such that:
116+ - If rho > 0, the outer eigenvalue lies in (d[n-1], infty), and f is increasing in this interval.
117+ - If rho < 0, the outer eigenvalue lies in (-infty, d[0]), and f is decreasing in this interval.
118+
119+ This function:
120+ 1. Determines the direction to search based on the sign of rho.
121+ 2. Finds an upper bound (or lower bound for rho < 0) where the secular function changes sign.
122+ 3. Uses the bisection method to find the root in the determined interval.
123+
124+ Parameters:
125+ f (callable): The secular function to find a root of.
126+ rho (float): Scalar rank-one update parameter.
127+ interval_end (float): Either d[0] or d[n-1], depending on the sign of rho.
128+ v (np.ndarray or scipy.sparse.spmatrix): Vector from the rank-one update; used to scale the search step size.
129+ tol (float, optional): Convergence tolerance. Default is 1e-12.
130+ max_iter (int, optional): Maximum number of bisection iterations. Default is 2000.
131+
132+ Returns:
133+ float: Approximation of the outer eigenvalue.
79134 """
80135 threshold = 1e-11
81136 update = np .linalg .norm (v )
@@ -97,10 +152,18 @@ def compute_outer_zero(f, rho, interval_end, v, tol=1e-12, max_iter=2000):
97152
98153def compute_psi_s (lambda_guess , rho , d , v_squared , i ):
99154 """
100- This function computes psi1, psi2 and their derivatives.
101- Corresponding to the interval (d[i], d[i+1]).
102- Unless I made some mistake, in case rho!=0, it should be sufficient to multiply all the sums by rho, which is what
103- is done here and seems to work in case rho > 0.
155+ Computes partial sums (psi1, psi2) and their derivatives for the secular function
156+ in the i-th interval (d[i], d[i+1]).
157+
158+ Parameters:
159+ lambda_guess (float): Evaluation point for the secular function.
160+ rho (float): Scalar rank-one update parameter.
161+ d (np.ndarray or scipy.sparse.spmatrix): 1D array of diagonal entries.
162+ v_squared (np.ndarray or scipy.sparse.spmatrix): Precomputed elementwise square of the update vector v.
163+ i (int): Index defining the interval (d[i], d[i+1]).
164+
165+ Returns:
166+ tuple: (psi1, psi1', psi2, psi2') — the partial secular sums and their derivatives.
104167 """
105168 denom1 = d [: i + 1 ] - lambda_guess
106169 denom2 = d [i + 1 :] - lambda_guess
@@ -113,7 +176,19 @@ def compute_psi_s(lambda_guess, rho, d, v_squared, i):
113176
114177def compute_inner_zero (rho , d , v , i , tol = 1e-12 , max_iter = 1000 ):
115178 """
116- This function (or some of its dependencies) surely contains a bug.
179+ Computes the i-th eigenvalue that lies in the interval (d[i], d[i+1]) for a
180+ rank-one modified diagonal matrix using the secular equation.
181+
182+ Parameters:
183+ rho (float): Rank-one update scalar.
184+ d (np.ndarray or scipy.sparse.spmatrix): 1D array of diagonal entries, sorted in ascending order.
185+ v (np.ndarray or scipy.sparse.spmatrix): 1D update vector.
186+ i (int): Index indicating the interval (d[i], d[i+1]) to find the zero in.
187+ tol (float, optional): Tolerance for root-finding. Default is 1e-12.
188+ max_iter (int, optional): Maximum iterations for bisection. Default is 1000.
189+
190+ Returns:
191+ float: The computed inner eigenvalue in the interval (d[i], d[i+1]).
117192 """
118193
119194 # fix the correct interval
@@ -152,6 +227,18 @@ def compute_inner_zero(rho, d, v, i, tol=1e-12, max_iter=1000):
152227
153228
154229def compute_eigenvalues (rho , d , v ):
230+ """
231+ Computes all eigenvalues of a rank-one modified diagonal matrix D + rho * v v^T
232+ using the secular equation method.
233+
234+ Parameters:
235+ rho (float): Rank-one update scalar.
236+ d (np.ndarray or scipy.sparse.spmatrix): 1D array of sorted diagonal entries of D.
237+ v (np.ndarray or scipy.sparse.spmatrix): Update vector v in the rank-one perturbation.
238+
239+ Returns:
240+ np.ndarray: Sorted array of all eigenvalues of the perturbed matrix.
241+ """
155242 f = return_secular_f (rho , d , v )
156243 eigenvalues = []
157244 iter_range = range (len (d ) - 1 )
0 commit comments