qdldl

import qdldl import numpy as np from scipy.sparse import csc_matrix # Construct a sparse symmetric indefinite matrix, typically a KKT matrix. # K = [ P A.T ] # [ A 0 ] # Example components: P_data = np.array([2.0, 1.0, 1.0, 3.0]) P_row_ind = np.array([0, 0, 1, 1]) P_col_ind = np.array([0, 1, 0, 1]) P = csc_matrix((P_data, (P_row_ind, P_col_ind)), shape=(2, 2)) A_constr_data = np.array([1.0, 1.0, 1.0]) A_constr_row_ind = np.array([0, 1, 1]) A_constr_col_ind = np.array([0, 0, 1]) A_constr = csc_matrix((A_constr_data, (A_constr_row_ind, A_constr_col_ind)), shape=(2, 2)) n_P = P.shape[0] n_A = A_constr.shape[0] n = n_P + n_A # Build the full KKT matrix KKT = csc_matrix( (n, n), dtype=P.dtype ) KKT[:n_P, :n_P] = P KKT[:n_P, n_P:] = A_constr.T KKT[n_P:, :n_P] = A_constr # The 'd' argument for qdldl.qdldl is a user-defined initial diagonal for pivoting. # A common choice is a vector of ones. d = np.ones(n) # Perform LDL factorization # L_factor is a csc_matrix (lower triangular factor), # D_diag is a numpy array (diagonal of D), # P_vec is a numpy array (permutation vector). L_factor, D_diag, P_vec = qdldl.qdldl(KKT, d) # Solve a system KKT @ x = b b = np.array([10.0, 20.0, 5.0, 8.0]) # Example right-hand side # Use the qdldl.solve function to get the solution x x = qdldl.solve(L_factor, D_diag, P_vec, b) print("Input KKT matrix (dense representation for display):\n", KKT.toarray()) print("Right-hand side b:\n", b) print("Solution x:\n", x) # Verify the solution print("KKT @ x:\n", KKT @ x) print("Error (KKT @ x - b):\n", KKT @ x - b) assert np.allclose(KKT @ x, b, atol=1e-9) print("Solution verified.")