Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_onenormest.py: 11%
198 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Sparse block 1-norm estimator.
2"""
4import numpy as np
5from scipy.sparse.linalg import aslinearoperator
8__all__ = ['onenormest']
11def onenormest(A, t=2, itmax=5, compute_v=False, compute_w=False):
12 """
13 Compute a lower bound of the 1-norm of a sparse matrix.
15 Parameters
16 ----------
17 A : ndarray or other linear operator
18 A linear operator that can be transposed and that can
19 produce matrix products.
20 t : int, optional
21 A positive parameter controlling the tradeoff between
22 accuracy versus time and memory usage.
23 Larger values take longer and use more memory
24 but give more accurate output.
25 itmax : int, optional
26 Use at most this many iterations.
27 compute_v : bool, optional
28 Request a norm-maximizing linear operator input vector if True.
29 compute_w : bool, optional
30 Request a norm-maximizing linear operator output vector if True.
32 Returns
33 -------
34 est : float
35 An underestimate of the 1-norm of the sparse matrix.
36 v : ndarray, optional
37 The vector such that ||Av||_1 == est*||v||_1.
38 It can be thought of as an input to the linear operator
39 that gives an output with particularly large norm.
40 w : ndarray, optional
41 The vector Av which has relatively large 1-norm.
42 It can be thought of as an output of the linear operator
43 that is relatively large in norm compared to the input.
45 Notes
46 -----
47 This is algorithm 2.4 of [1].
49 In [2] it is described as follows.
50 "This algorithm typically requires the evaluation of
51 about 4t matrix-vector products and almost invariably
52 produces a norm estimate (which is, in fact, a lower
53 bound on the norm) correct to within a factor 3."
55 .. versionadded:: 0.13.0
57 References
58 ----------
59 .. [1] Nicholas J. Higham and Francoise Tisseur (2000),
60 "A Block Algorithm for Matrix 1-Norm Estimation,
61 with an Application to 1-Norm Pseudospectra."
62 SIAM J. Matrix Anal. Appl. Vol. 21, No. 4, pp. 1185-1201.
64 .. [2] Awad H. Al-Mohy and Nicholas J. Higham (2009),
65 "A new scaling and squaring algorithm for the matrix exponential."
66 SIAM J. Matrix Anal. Appl. Vol. 31, No. 3, pp. 970-989.
68 Examples
69 --------
70 >>> import numpy as np
71 >>> from scipy.sparse import csc_matrix
72 >>> from scipy.sparse.linalg import onenormest
73 >>> A = csc_matrix([[1., 0., 0.], [5., 8., 2.], [0., -1., 0.]], dtype=float)
74 >>> A.toarray()
75 array([[ 1., 0., 0.],
76 [ 5., 8., 2.],
77 [ 0., -1., 0.]])
78 >>> onenormest(A)
79 9.0
80 >>> np.linalg.norm(A.toarray(), ord=1)
81 9.0
82 """
84 # Check the input.
85 A = aslinearoperator(A)
86 if A.shape[0] != A.shape[1]:
87 raise ValueError('expected the operator to act like a square matrix')
89 # If the operator size is small compared to t,
90 # then it is easier to compute the exact norm.
91 # Otherwise estimate the norm.
92 n = A.shape[1]
93 if t >= n:
94 A_explicit = np.asarray(aslinearoperator(A).matmat(np.identity(n)))
95 if A_explicit.shape != (n, n):
96 raise Exception('internal error: ',
97 'unexpected shape ' + str(A_explicit.shape))
98 col_abs_sums = abs(A_explicit).sum(axis=0)
99 if col_abs_sums.shape != (n, ):
100 raise Exception('internal error: ',
101 'unexpected shape ' + str(col_abs_sums.shape))
102 argmax_j = np.argmax(col_abs_sums)
103 v = elementary_vector(n, argmax_j)
104 w = A_explicit[:, argmax_j]
105 est = col_abs_sums[argmax_j]
106 else:
107 est, v, w, nmults, nresamples = _onenormest_core(A, A.H, t, itmax)
109 # Report the norm estimate along with some certificates of the estimate.
110 if compute_v or compute_w:
111 result = (est,)
112 if compute_v:
113 result += (v,)
114 if compute_w:
115 result += (w,)
116 return result
117 else:
118 return est
121def _blocked_elementwise(func):
122 """
123 Decorator for an elementwise function, to apply it blockwise along
124 first dimension, to avoid excessive memory usage in temporaries.
125 """
126 block_size = 2**20
128 def wrapper(x):
129 if x.shape[0] < block_size:
130 return func(x)
131 else:
132 y0 = func(x[:block_size])
133 y = np.zeros((x.shape[0],) + y0.shape[1:], dtype=y0.dtype)
134 y[:block_size] = y0
135 del y0
136 for j in range(block_size, x.shape[0], block_size):
137 y[j:j+block_size] = func(x[j:j+block_size])
138 return y
139 return wrapper
142@_blocked_elementwise
143def sign_round_up(X):
144 """
145 This should do the right thing for both real and complex matrices.
147 From Higham and Tisseur:
148 "Everything in this section remains valid for complex matrices
149 provided that sign(A) is redefined as the matrix (aij / |aij|)
150 (and sign(0) = 1) transposes are replaced by conjugate transposes."
152 """
153 Y = X.copy()
154 Y[Y == 0] = 1
155 Y /= np.abs(Y)
156 return Y
159@_blocked_elementwise
160def _max_abs_axis1(X):
161 return np.max(np.abs(X), axis=1)
164def _sum_abs_axis0(X):
165 block_size = 2**20
166 r = None
167 for j in range(0, X.shape[0], block_size):
168 y = np.sum(np.abs(X[j:j+block_size]), axis=0)
169 if r is None:
170 r = y
171 else:
172 r += y
173 return r
176def elementary_vector(n, i):
177 v = np.zeros(n, dtype=float)
178 v[i] = 1
179 return v
182def vectors_are_parallel(v, w):
183 # Columns are considered parallel when they are equal or negative.
184 # Entries are required to be in {-1, 1},
185 # which guarantees that the magnitudes of the vectors are identical.
186 if v.ndim != 1 or v.shape != w.shape:
187 raise ValueError('expected conformant vectors with entries in {-1,1}')
188 n = v.shape[0]
189 return np.dot(v, w) == n
192def every_col_of_X_is_parallel_to_a_col_of_Y(X, Y):
193 for v in X.T:
194 if not any(vectors_are_parallel(v, w) for w in Y.T):
195 return False
196 return True
199def column_needs_resampling(i, X, Y=None):
200 # column i of X needs resampling if either
201 # it is parallel to a previous column of X or
202 # it is parallel to a column of Y
203 n, t = X.shape
204 v = X[:, i]
205 if any(vectors_are_parallel(v, X[:, j]) for j in range(i)):
206 return True
207 if Y is not None:
208 if any(vectors_are_parallel(v, w) for w in Y.T):
209 return True
210 return False
213def resample_column(i, X):
214 X[:, i] = np.random.randint(0, 2, size=X.shape[0])*2 - 1
217def less_than_or_close(a, b):
218 return np.allclose(a, b) or (a < b)
221def _algorithm_2_2(A, AT, t):
222 """
223 This is Algorithm 2.2.
225 Parameters
226 ----------
227 A : ndarray or other linear operator
228 A linear operator that can produce matrix products.
229 AT : ndarray or other linear operator
230 The transpose of A.
231 t : int, optional
232 A positive parameter controlling the tradeoff between
233 accuracy versus time and memory usage.
235 Returns
236 -------
237 g : sequence
238 A non-negative decreasing vector
239 such that g[j] is a lower bound for the 1-norm
240 of the column of A of jth largest 1-norm.
241 The first entry of this vector is therefore a lower bound
242 on the 1-norm of the linear operator A.
243 This sequence has length t.
244 ind : sequence
245 The ith entry of ind is the index of the column A whose 1-norm
246 is given by g[i].
247 This sequence of indices has length t, and its entries are
248 chosen from range(n), possibly with repetition,
249 where n is the order of the operator A.
251 Notes
252 -----
253 This algorithm is mainly for testing.
254 It uses the 'ind' array in a way that is similar to
255 its usage in algorithm 2.4. This algorithm 2.2 may be easier to test,
256 so it gives a chance of uncovering bugs related to indexing
257 which could have propagated less noticeably to algorithm 2.4.
259 """
260 A_linear_operator = aslinearoperator(A)
261 AT_linear_operator = aslinearoperator(AT)
262 n = A_linear_operator.shape[0]
264 # Initialize the X block with columns of unit 1-norm.
265 X = np.ones((n, t))
266 if t > 1:
267 X[:, 1:] = np.random.randint(0, 2, size=(n, t-1))*2 - 1
268 X /= float(n)
270 # Iteratively improve the lower bounds.
271 # Track extra things, to assert invariants for debugging.
272 g_prev = None
273 h_prev = None
274 k = 1
275 ind = range(t)
276 while True:
277 Y = np.asarray(A_linear_operator.matmat(X))
278 g = _sum_abs_axis0(Y)
279 best_j = np.argmax(g)
280 g.sort()
281 g = g[::-1]
282 S = sign_round_up(Y)
283 Z = np.asarray(AT_linear_operator.matmat(S))
284 h = _max_abs_axis1(Z)
286 # If this algorithm runs for fewer than two iterations,
287 # then its return values do not have the properties indicated
288 # in the description of the algorithm.
289 # In particular, the entries of g are not 1-norms of any
290 # column of A until the second iteration.
291 # Therefore we will require the algorithm to run for at least
292 # two iterations, even though this requirement is not stated
293 # in the description of the algorithm.
294 if k >= 2:
295 if less_than_or_close(max(h), np.dot(Z[:, best_j], X[:, best_j])):
296 break
297 ind = np.argsort(h)[::-1][:t]
298 h = h[ind]
299 for j in range(t):
300 X[:, j] = elementary_vector(n, ind[j])
302 # Check invariant (2.2).
303 if k >= 2:
304 if not less_than_or_close(g_prev[0], h_prev[0]):
305 raise Exception('invariant (2.2) is violated')
306 if not less_than_or_close(h_prev[0], g[0]):
307 raise Exception('invariant (2.2) is violated')
309 # Check invariant (2.3).
310 if k >= 3:
311 for j in range(t):
312 if not less_than_or_close(g[j], g_prev[j]):
313 raise Exception('invariant (2.3) is violated')
315 # Update for the next iteration.
316 g_prev = g
317 h_prev = h
318 k += 1
320 # Return the lower bounds and the corresponding column indices.
321 return g, ind
324def _onenormest_core(A, AT, t, itmax):
325 """
326 Compute a lower bound of the 1-norm of a sparse matrix.
328 Parameters
329 ----------
330 A : ndarray or other linear operator
331 A linear operator that can produce matrix products.
332 AT : ndarray or other linear operator
333 The transpose of A.
334 t : int, optional
335 A positive parameter controlling the tradeoff between
336 accuracy versus time and memory usage.
337 itmax : int, optional
338 Use at most this many iterations.
340 Returns
341 -------
342 est : float
343 An underestimate of the 1-norm of the sparse matrix.
344 v : ndarray, optional
345 The vector such that ||Av||_1 == est*||v||_1.
346 It can be thought of as an input to the linear operator
347 that gives an output with particularly large norm.
348 w : ndarray, optional
349 The vector Av which has relatively large 1-norm.
350 It can be thought of as an output of the linear operator
351 that is relatively large in norm compared to the input.
352 nmults : int, optional
353 The number of matrix products that were computed.
354 nresamples : int, optional
355 The number of times a parallel column was observed,
356 necessitating a re-randomization of the column.
358 Notes
359 -----
360 This is algorithm 2.4.
362 """
363 # This function is a more or less direct translation
364 # of Algorithm 2.4 from the Higham and Tisseur (2000) paper.
365 A_linear_operator = aslinearoperator(A)
366 AT_linear_operator = aslinearoperator(AT)
367 if itmax < 2:
368 raise ValueError('at least two iterations are required')
369 if t < 1:
370 raise ValueError('at least one column is required')
371 n = A.shape[0]
372 if t >= n:
373 raise ValueError('t should be smaller than the order of A')
374 # Track the number of big*small matrix multiplications
375 # and the number of resamplings.
376 nmults = 0
377 nresamples = 0
378 # "We now explain our choice of starting matrix. We take the first
379 # column of X to be the vector of 1s [...] This has the advantage that
380 # for a matrix with nonnegative elements the algorithm converges
381 # with an exact estimate on the second iteration, and such matrices
382 # arise in applications [...]"
383 X = np.ones((n, t), dtype=float)
384 # "The remaining columns are chosen as rand{-1,1},
385 # with a check for and correction of parallel columns,
386 # exactly as for S in the body of the algorithm."
387 if t > 1:
388 for i in range(1, t):
389 # These are technically initial samples, not resamples,
390 # so the resampling count is not incremented.
391 resample_column(i, X)
392 for i in range(t):
393 while column_needs_resampling(i, X):
394 resample_column(i, X)
395 nresamples += 1
396 # "Choose starting matrix X with columns of unit 1-norm."
397 X /= float(n)
398 # "indices of used unit vectors e_j"
399 ind_hist = np.zeros(0, dtype=np.intp)
400 est_old = 0
401 S = np.zeros((n, t), dtype=float)
402 k = 1
403 ind = None
404 while True:
405 Y = np.asarray(A_linear_operator.matmat(X))
406 nmults += 1
407 mags = _sum_abs_axis0(Y)
408 est = np.max(mags)
409 best_j = np.argmax(mags)
410 if est > est_old or k == 2:
411 if k >= 2:
412 ind_best = ind[best_j]
413 w = Y[:, best_j]
414 # (1)
415 if k >= 2 and est <= est_old:
416 est = est_old
417 break
418 est_old = est
419 S_old = S
420 if k > itmax:
421 break
422 S = sign_round_up(Y)
423 del Y
424 # (2)
425 if every_col_of_X_is_parallel_to_a_col_of_Y(S, S_old):
426 break
427 if t > 1:
428 # "Ensure that no column of S is parallel to another column of S
429 # or to a column of S_old by replacing columns of S by rand{-1,1}."
430 for i in range(t):
431 while column_needs_resampling(i, S, S_old):
432 resample_column(i, S)
433 nresamples += 1
434 del S_old
435 # (3)
436 Z = np.asarray(AT_linear_operator.matmat(S))
437 nmults += 1
438 h = _max_abs_axis1(Z)
439 del Z
440 # (4)
441 if k >= 2 and max(h) == h[ind_best]:
442 break
443 # "Sort h so that h_first >= ... >= h_last
444 # and re-order ind correspondingly."
445 #
446 # Later on, we will need at most t+len(ind_hist) largest
447 # entries, so drop the rest
448 ind = np.argsort(h)[::-1][:t+len(ind_hist)].copy()
449 del h
450 if t > 1:
451 # (5)
452 # Break if the most promising t vectors have been visited already.
453 if np.in1d(ind[:t], ind_hist).all():
454 break
455 # Put the most promising unvisited vectors at the front of the list
456 # and put the visited vectors at the end of the list.
457 # Preserve the order of the indices induced by the ordering of h.
458 seen = np.in1d(ind, ind_hist)
459 ind = np.concatenate((ind[~seen], ind[seen]))
460 for j in range(t):
461 X[:, j] = elementary_vector(n, ind[j])
463 new_ind = ind[:t][~np.in1d(ind[:t], ind_hist)]
464 ind_hist = np.concatenate((ind_hist, new_ind))
465 k += 1
466 v = elementary_vector(n, ind_best)
467 return est, v, w, nmults, nresamples