|
17 | 17 |
|
18 | 18 | import numpy as np |
19 | 19 | import scipy.sparse as sp |
20 | | -from scipy import linalg as LA |
21 | 20 |
|
22 | 21 | from cvxpy.atoms.affine.wraps import psd_wrap |
23 | 22 | from cvxpy.atoms.atom import Atom |
24 | 23 | from cvxpy.expressions.constants.parameter import is_param_affine, is_param_free |
25 | 24 | from cvxpy.expressions.expression import Expression |
26 | 25 | from cvxpy.interface.matrix_utilities import is_sparse |
27 | 26 | from cvxpy.utilities import scopes |
28 | | -from cvxpy.utilities.linalg import sparse_cholesky |
| 27 | +from cvxpy.utilities.linalg import dense_ldl_decomp, sparse_cholesky |
29 | 28 | from cvxpy.utilities.warn import warn |
30 | 29 |
|
31 | 30 |
|
@@ -264,27 +263,7 @@ def decomp_quad(P, cond=None, rcond=None, lower=True, check_finite: bool = True) |
264 | 263 | return 1.0, np.empty((0, 0)), L[p, :] |
265 | 264 | except ValueError: |
266 | 265 | P = P.toarray() # make dense (needs to happen for ldl). |
267 | | - lu, d, _perm = LA.ldl(P, lower=lower, check_finite=check_finite) |
268 | | - |
269 | | - # Extract effective diagonal values from D, handling any 2x2 blocks. |
270 | | - # For PSD/NSD matrices D is diagonal (no 2x2 blocks). For indefinite |
271 | | - # matrices, Bunch-Kaufman pivoting may introduce 2x2 blocks which we |
272 | | - # resolve by batching all 2x2 blocks into a single eigh call. |
273 | | - sub_diag = np.diag(d, -1) |
274 | | - block_starts = np.nonzero(sub_diag)[0] |
275 | | - diag_vals = np.real(np.diag(d)).copy() |
276 | | - if len(block_starts) > 0: |
277 | | - bs = block_starts |
278 | | - idx = bs[:, None] + np.arange(2)[None, :] # (k, 2) |
279 | | - blocks = d[idx[:, :, None], idx[:, None, :]] # (k, 2, 2) |
280 | | - eigvals, eigvecs = np.linalg.eigh(blocks) # (k, 2), (k, 2, 2) |
281 | | - diag_vals[bs] = eigvals[:, 0] |
282 | | - diag_vals[bs + 1] = eigvals[:, 1] |
283 | | - # Apply eigenvector rotations to lu columns |
284 | | - lu_i = lu[:, bs].copy() |
285 | | - lu_ip1 = lu[:, bs + 1].copy() |
286 | | - lu[:, bs] = lu_i * eigvecs[:, 0, 0] + lu_ip1 * eigvecs[:, 1, 0] |
287 | | - lu[:, bs + 1] = lu_i * eigvecs[:, 0, 1] + lu_ip1 * eigvecs[:, 1, 1] |
| 266 | + diag_vals, lu = dense_ldl_decomp(P, lower=lower, check_finite=check_finite) |
288 | 267 |
|
289 | 268 | if rcond is not None: |
290 | 269 | cond = rcond |
|
0 commit comments