Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_remove_redundancy.py: 8%
173 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2Routines for removing redundant (linearly dependent) equations from linear
3programming equality constraints.
4"""
5# Author: Matt Haberland
7import numpy as np
8from scipy.linalg import svd
9from scipy.linalg.interpolative import interp_decomp
10import scipy
11from scipy.linalg.blas import dtrsm
14def _row_count(A):
15 """
16 Counts the number of nonzeros in each row of input array A.
17 Nonzeros are defined as any element with absolute value greater than
18 tol = 1e-13. This value should probably be an input to the function.
20 Parameters
21 ----------
22 A : 2-D array
23 An array representing a matrix
25 Returns
26 -------
27 rowcount : 1-D array
28 Number of nonzeros in each row of A
30 """
31 tol = 1e-13
32 return np.array((abs(A) > tol).sum(axis=1)).flatten()
35def _get_densest(A, eligibleRows):
36 """
37 Returns the index of the densest row of A. Ignores rows that are not
38 eligible for consideration.
40 Parameters
41 ----------
42 A : 2-D array
43 An array representing a matrix
44 eligibleRows : 1-D logical array
45 Values indicate whether the corresponding row of A is eligible
46 to be considered
48 Returns
49 -------
50 i_densest : int
51 Index of the densest row in A eligible for consideration
53 """
54 rowCounts = _row_count(A)
55 return np.argmax(rowCounts * eligibleRows)
58def _remove_zero_rows(A, b):
59 """
60 Eliminates trivial equations from system of equations defined by Ax = b
61 and identifies trivial infeasibilities
63 Parameters
64 ----------
65 A : 2-D array
66 An array representing the left-hand side of a system of equations
67 b : 1-D array
68 An array representing the right-hand side of a system of equations
70 Returns
71 -------
72 A : 2-D array
73 An array representing the left-hand side of a system of equations
74 b : 1-D array
75 An array representing the right-hand side of a system of equations
76 status: int
77 An integer indicating the status of the removal operation
78 0: No infeasibility identified
79 2: Trivially infeasible
80 message : str
81 A string descriptor of the exit status of the optimization.
83 """
84 status = 0
85 message = ""
86 i_zero = _row_count(A) == 0
87 A = A[np.logical_not(i_zero), :]
88 if not np.allclose(b[i_zero], 0):
89 status = 2
90 message = "There is a zero row in A_eq with a nonzero corresponding " \
91 "entry in b_eq. The problem is infeasible."
92 b = b[np.logical_not(i_zero)]
93 return A, b, status, message
96def bg_update_dense(plu, perm_r, v, j):
97 LU, p = plu
99 vperm = v[perm_r]
100 u = dtrsm(1, LU, vperm, lower=1, diag=1)
101 LU[:j+1, j] = u[:j+1]
102 l = u[j+1:]
103 piv = LU[j, j]
104 LU[j+1:, j] += (l/piv)
105 return LU, p
108def _remove_redundancy_pivot_dense(A, rhs, true_rank=None):
109 """
110 Eliminates redundant equations from system of equations defined by Ax = b
111 and identifies infeasibilities.
113 Parameters
114 ----------
115 A : 2-D sparse matrix
116 An matrix representing the left-hand side of a system of equations
117 rhs : 1-D array
118 An array representing the right-hand side of a system of equations
120 Returns
121 -------
122 A : 2-D sparse matrix
123 A matrix representing the left-hand side of a system of equations
124 rhs : 1-D array
125 An array representing the right-hand side of a system of equations
126 status: int
127 An integer indicating the status of the system
128 0: No infeasibility identified
129 2: Trivially infeasible
130 message : str
131 A string descriptor of the exit status of the optimization.
133 References
134 ----------
135 .. [2] Andersen, Erling D. "Finding all linearly dependent rows in
136 large-scale linear programming." Optimization Methods and Software
137 6.3 (1995): 219-227.
139 """
140 tolapiv = 1e-8
141 tolprimal = 1e-8
142 status = 0
143 message = ""
144 inconsistent = ("There is a linear combination of rows of A_eq that "
145 "results in zero, suggesting a redundant constraint. "
146 "However the same linear combination of b_eq is "
147 "nonzero, suggesting that the constraints conflict "
148 "and the problem is infeasible.")
149 A, rhs, status, message = _remove_zero_rows(A, rhs)
151 if status != 0:
152 return A, rhs, status, message
154 m, n = A.shape
156 v = list(range(m)) # Artificial column indices.
157 b = list(v) # Basis column indices.
158 # This is better as a list than a set because column order of basis matrix
159 # needs to be consistent.
160 d = [] # Indices of dependent rows
161 perm_r = None
163 A_orig = A
164 A = np.zeros((m, m + n), order='F')
165 np.fill_diagonal(A, 1)
166 A[:, m:] = A_orig
167 e = np.zeros(m)
169 js_candidates = np.arange(m, m+n, dtype=int) # candidate columns for basis
170 # manual masking was faster than masked array
171 js_mask = np.ones(js_candidates.shape, dtype=bool)
173 # Implements basic algorithm from [2]
174 # Uses some of the suggested improvements (removing zero rows and
175 # Bartels-Golub update idea).
176 # Removing column singletons would be easy, but it is not as important
177 # because the procedure is performed only on the equality constraint
178 # matrix from the original problem - not on the canonical form matrix,
179 # which would have many more column singletons due to slack variables
180 # from the inequality constraints.
181 # The thoughts on "crashing" the initial basis are only really useful if
182 # the matrix is sparse.
184 lu = np.eye(m, order='F'), np.arange(m) # initial LU is trivial
185 perm_r = lu[1]
186 for i in v:
188 e[i] = 1
189 if i > 0:
190 e[i-1] = 0
192 try: # fails for i==0 and any time it gets ill-conditioned
193 j = b[i-1]
194 lu = bg_update_dense(lu, perm_r, A[:, j], i-1)
195 except Exception:
196 lu = scipy.linalg.lu_factor(A[:, b])
197 LU, p = lu
198 perm_r = list(range(m))
199 for i1, i2 in enumerate(p):
200 perm_r[i1], perm_r[i2] = perm_r[i2], perm_r[i1]
202 pi = scipy.linalg.lu_solve(lu, e, trans=1)
204 js = js_candidates[js_mask]
205 batch = 50
207 # This is a tiny bit faster than looping over columns indivually,
208 # like for j in js: if abs(A[:,j].transpose().dot(pi)) > tolapiv:
209 for j_index in range(0, len(js), batch):
210 j_indices = js[j_index: min(j_index+batch, len(js))]
212 c = abs(A[:, j_indices].transpose().dot(pi))
213 if (c > tolapiv).any():
214 j = js[j_index + np.argmax(c)] # very independent column
215 b[i] = j
216 js_mask[j-m] = False
217 break
218 else:
219 bibar = pi.T.dot(rhs.reshape(-1, 1))
220 bnorm = np.linalg.norm(rhs)
221 if abs(bibar)/(1+bnorm) > tolprimal: # inconsistent
222 status = 2
223 message = inconsistent
224 return A_orig, rhs, status, message
225 else: # dependent
226 d.append(i)
227 if true_rank is not None and len(d) == m - true_rank:
228 break # found all redundancies
230 keep = set(range(m))
231 keep = list(keep - set(d))
232 return A_orig[keep, :], rhs[keep], status, message
235def _remove_redundancy_pivot_sparse(A, rhs):
236 """
237 Eliminates redundant equations from system of equations defined by Ax = b
238 and identifies infeasibilities.
240 Parameters
241 ----------
242 A : 2-D sparse matrix
243 An matrix representing the left-hand side of a system of equations
244 rhs : 1-D array
245 An array representing the right-hand side of a system of equations
247 Returns
248 -------
249 A : 2-D sparse matrix
250 A matrix representing the left-hand side of a system of equations
251 rhs : 1-D array
252 An array representing the right-hand side of a system of equations
253 status: int
254 An integer indicating the status of the system
255 0: No infeasibility identified
256 2: Trivially infeasible
257 message : str
258 A string descriptor of the exit status of the optimization.
260 References
261 ----------
262 .. [2] Andersen, Erling D. "Finding all linearly dependent rows in
263 large-scale linear programming." Optimization Methods and Software
264 6.3 (1995): 219-227.
266 """
268 tolapiv = 1e-8
269 tolprimal = 1e-8
270 status = 0
271 message = ""
272 inconsistent = ("There is a linear combination of rows of A_eq that "
273 "results in zero, suggesting a redundant constraint. "
274 "However the same linear combination of b_eq is "
275 "nonzero, suggesting that the constraints conflict "
276 "and the problem is infeasible.")
277 A, rhs, status, message = _remove_zero_rows(A, rhs)
279 if status != 0:
280 return A, rhs, status, message
282 m, n = A.shape
284 v = list(range(m)) # Artificial column indices.
285 b = list(v) # Basis column indices.
286 # This is better as a list than a set because column order of basis matrix
287 # needs to be consistent.
288 k = set(range(m, m+n)) # Structural column indices.
289 d = [] # Indices of dependent rows
291 A_orig = A
292 A = scipy.sparse.hstack((scipy.sparse.eye(m), A)).tocsc()
293 e = np.zeros(m)
295 # Implements basic algorithm from [2]
296 # Uses only one of the suggested improvements (removing zero rows).
297 # Removing column singletons would be easy, but it is not as important
298 # because the procedure is performed only on the equality constraint
299 # matrix from the original problem - not on the canonical form matrix,
300 # which would have many more column singletons due to slack variables
301 # from the inequality constraints.
302 # The thoughts on "crashing" the initial basis sound useful, but the
303 # description of the procedure seems to assume a lot of familiarity with
304 # the subject; it is not very explicit. I already went through enough
305 # trouble getting the basic algorithm working, so I was not interested in
306 # trying to decipher this, too. (Overall, the paper is fraught with
307 # mistakes and ambiguities - which is strange, because the rest of
308 # Andersen's papers are quite good.)
309 # I tried and tried and tried to improve performance using the
310 # Bartels-Golub update. It works, but it's only practical if the LU
311 # factorization can be specialized as described, and that is not possible
312 # until the SciPy SuperLU interface permits control over column
313 # permutation - see issue #7700.
315 for i in v:
316 B = A[:, b]
318 e[i] = 1
319 if i > 0:
320 e[i-1] = 0
322 pi = scipy.sparse.linalg.spsolve(B.transpose(), e).reshape(-1, 1)
324 js = list(k-set(b)) # not efficient, but this is not the time sink...
326 # Due to overhead, it tends to be faster (for problems tested) to
327 # compute the full matrix-vector product rather than individual
328 # vector-vector products (with the chance of terminating as soon
329 # as any are nonzero). For very large matrices, it might be worth
330 # it to compute, say, 100 or 1000 at a time and stop when a nonzero
331 # is found.
333 c = (np.abs(A[:, js].transpose().dot(pi)) > tolapiv).nonzero()[0]
334 if len(c) > 0: # independent
335 j = js[c[0]]
336 # in a previous commit, the previous line was changed to choose
337 # index j corresponding with the maximum dot product.
338 # While this avoided issues with almost
339 # singular matrices, it slowed the routine in most NETLIB tests.
340 # I think this is because these columns were denser than the
341 # first column with nonzero dot product (c[0]).
342 # It would be nice to have a heuristic that balances sparsity with
343 # high dot product, but I don't think it's worth the time to
344 # develop one right now. Bartels-Golub update is a much higher
345 # priority.
346 b[i] = j # replace artificial column
347 else:
348 bibar = pi.T.dot(rhs.reshape(-1, 1))
349 bnorm = np.linalg.norm(rhs)
350 if abs(bibar)/(1 + bnorm) > tolprimal:
351 status = 2
352 message = inconsistent
353 return A_orig, rhs, status, message
354 else: # dependent
355 d.append(i)
357 keep = set(range(m))
358 keep = list(keep - set(d))
359 return A_orig[keep, :], rhs[keep], status, message
362def _remove_redundancy_svd(A, b):
363 """
364 Eliminates redundant equations from system of equations defined by Ax = b
365 and identifies infeasibilities.
367 Parameters
368 ----------
369 A : 2-D array
370 An array representing the left-hand side of a system of equations
371 b : 1-D array
372 An array representing the right-hand side of a system of equations
374 Returns
375 -------
376 A : 2-D array
377 An array representing the left-hand side of a system of equations
378 b : 1-D array
379 An array representing the right-hand side of a system of equations
380 status: int
381 An integer indicating the status of the system
382 0: No infeasibility identified
383 2: Trivially infeasible
384 message : str
385 A string descriptor of the exit status of the optimization.
387 References
388 ----------
389 .. [2] Andersen, Erling D. "Finding all linearly dependent rows in
390 large-scale linear programming." Optimization Methods and Software
391 6.3 (1995): 219-227.
393 """
395 A, b, status, message = _remove_zero_rows(A, b)
397 if status != 0:
398 return A, b, status, message
400 U, s, Vh = svd(A)
401 eps = np.finfo(float).eps
402 tol = s.max() * max(A.shape) * eps
404 m, n = A.shape
405 s_min = s[-1] if m <= n else 0
407 # this algorithm is faster than that of [2] when the nullspace is small
408 # but it could probably be improvement by randomized algorithms and with
409 # a sparse implementation.
410 # it relies on repeated singular value decomposition to find linearly
411 # dependent rows (as identified by columns of U that correspond with zero
412 # singular values). Unfortunately, only one row can be removed per
413 # decomposition (I tried otherwise; doing so can cause problems.)
414 # It would be nice if we could do truncated SVD like sp.sparse.linalg.svds
415 # but that function is unreliable at finding singular values near zero.
416 # Finding max eigenvalue L of A A^T, then largest eigenvalue (and
417 # associated eigenvector) of -A A^T + L I (I is identity) via power
418 # iteration would also work in theory, but is only efficient if the
419 # smallest nonzero eigenvalue of A A^T is close to the largest nonzero
420 # eigenvalue.
422 while abs(s_min) < tol:
423 v = U[:, -1] # TODO: return these so user can eliminate from problem?
424 # rows need to be represented in significant amount
425 eligibleRows = np.abs(v) > tol * 10e6
426 if not np.any(eligibleRows) or np.any(np.abs(v.dot(A)) > tol):
427 status = 4
428 message = ("Due to numerical issues, redundant equality "
429 "constraints could not be removed automatically. "
430 "Try providing your constraint matrices as sparse "
431 "matrices to activate sparse presolve, try turning "
432 "off redundancy removal, or try turning off presolve "
433 "altogether.")
434 break
435 if np.any(np.abs(v.dot(b)) > tol * 100): # factor of 100 to fix 10038 and 10349
436 status = 2
437 message = ("There is a linear combination of rows of A_eq that "
438 "results in zero, suggesting a redundant constraint. "
439 "However the same linear combination of b_eq is "
440 "nonzero, suggesting that the constraints conflict "
441 "and the problem is infeasible.")
442 break
444 i_remove = _get_densest(A, eligibleRows)
445 A = np.delete(A, i_remove, axis=0)
446 b = np.delete(b, i_remove)
447 U, s, Vh = svd(A)
448 m, n = A.shape
449 s_min = s[-1] if m <= n else 0
451 return A, b, status, message
454def _remove_redundancy_id(A, rhs, rank=None, randomized=True):
455 """Eliminates redundant equations from a system of equations.
457 Eliminates redundant equations from system of equations defined by Ax = b
458 and identifies infeasibilities.
460 Parameters
461 ----------
462 A : 2-D array
463 An array representing the left-hand side of a system of equations
464 rhs : 1-D array
465 An array representing the right-hand side of a system of equations
466 rank : int, optional
467 The rank of A
468 randomized: bool, optional
469 True for randomized interpolative decomposition
471 Returns
472 -------
473 A : 2-D array
474 An array representing the left-hand side of a system of equations
475 rhs : 1-D array
476 An array representing the right-hand side of a system of equations
477 status: int
478 An integer indicating the status of the system
479 0: No infeasibility identified
480 2: Trivially infeasible
481 message : str
482 A string descriptor of the exit status of the optimization.
484 """
486 status = 0
487 message = ""
488 inconsistent = ("There is a linear combination of rows of A_eq that "
489 "results in zero, suggesting a redundant constraint. "
490 "However the same linear combination of b_eq is "
491 "nonzero, suggesting that the constraints conflict "
492 "and the problem is infeasible.")
494 A, rhs, status, message = _remove_zero_rows(A, rhs)
496 if status != 0:
497 return A, rhs, status, message
499 m, n = A.shape
501 k = rank
502 if rank is None:
503 k = np.linalg.matrix_rank(A)
505 idx, proj = interp_decomp(A.T, k, rand=randomized)
507 # first k entries in idx are indices of the independent rows
508 # remaining entries are the indices of the m-k dependent rows
509 # proj provides a linear combinations of rows of A2 that form the
510 # remaining m-k (dependent) rows. The same linear combination of entries
511 # in rhs2 must give the remaining m-k entries. If not, the system is
512 # inconsistent, and the problem is infeasible.
513 if not np.allclose(rhs[idx[:k]] @ proj, rhs[idx[k:]]):
514 status = 2
515 message = inconsistent
517 # sort indices because the other redundancy removal routines leave rows
518 # in original order and tests were written with that in mind
519 idx = sorted(idx[:k])
520 A2 = A[idx, :]
521 rhs2 = rhs[idx]
522 return A2, rhs2, status, message