From 89b379d05f51b926573031c3f5a3b9c06769bf87 Mon Sep 17 00:00:00 2001 From: darnstrom Date: Fri, 14 Aug 2026 21:06:05 +0200 Subject: [PATCH] Detect symmetric AVIs --- include/types.h | 2 + include/utils.h | 2 +- interfaces/daqp-python/test/example_test.py | 41 +++++++++++++++++++++ src/api.c | 8 +++- src/eq_elim.c | 3 +- src/utils.c | 34 ++++++++++++----- 6 files changed, 77 insertions(+), 13 deletions(-) diff --git a/include/types.h b/include/types.h index bd673d1..902e1d1 100644 --- a/include/types.h +++ b/include/types.h @@ -99,6 +99,8 @@ typedef struct{ }DAQPBnB; typedef struct{ + int is_symmetric; + c_float* Hsym; c_float* Hs_rho; c_float* H_rho; diff --git a/include/utils.h b/include/utils.h index 6d1d885..1ede787 100644 --- a/include/utils.h +++ b/include/utils.h @@ -19,7 +19,7 @@ void daqp_normalize_Rinv(DAQPWorkspace *work); int daqp_normalize_M(DAQPWorkspace *work); int daqp_check_unconstrained(DAQPWorkspace* work, const int mask); -int daqp_update_avi(DAQPAVI *avi, DAQPProblem *problem); +int daqp_update_avi(DAQPAVI *avi, DAQPProblem *problem, c_float zero_tol); int daqp_lu(c_float* A, int* P, int n); void daqp_lu_solve(c_float* LU, int* P, c_float* b, c_float* x, int n); diff --git a/interfaces/daqp-python/test/example_test.py b/interfaces/daqp-python/test/example_test.py index 6ed6f0d..308918c 100644 --- a/interfaces/daqp-python/test/example_test.py +++ b/interfaces/daqp-python/test/example_test.py @@ -25,6 +25,47 @@ def test_python_demo(self): xstar, fval, exitflag, lam = daqp.solve(H, f, A, bupper, blower, sense) self.assertEqual(exitflag, 1) + def test_symmetric_avi_uses_qp_fast_path(self): + """An ill-conditioned symmetric AVI is solved as its equivalent QP.""" + Q = np.array([[1.0, 1.0], [-1.0, 1.0]], dtype=c_double) / np.sqrt(2.0) + H = Q @ np.diag([1.0e-6, 1.0e6]) @ Q.T + f = np.array([1.0, 2.0], dtype=c_double) + A = np.eye(2, dtype=c_double) + bupper = np.ones(2, dtype=c_double) + blower = -np.ones(2, dtype=c_double) + sense = np.zeros(2, dtype=c_int) + + x_qp, _, flag_qp, _ = daqp.solve( + H, f, A, bupper, blower, sense, iter_limit=10) + x_avi, _, flag_avi, _ = daqp.solve( + H, f, A, bupper, blower, sense, is_avi=True, iter_limit=10) + + self.assertEqual(flag_qp, 1) + self.assertEqual(flag_avi, 1) + np.testing.assert_allclose(x_avi, x_qp, atol=1.0e-8) + + def test_avi_model_dispatches_on_symmetry(self): + """Persistent AVI models select the QP or DR path during setup.""" + H_symmetric = np.array([[2.0, 0.5], [0.5, 1.0]], dtype=c_double) + H_asymmetric = np.array([[2.0, 1.0], [0.0, 1.0]], dtype=c_double) + f = np.array([1.0, 2.0], dtype=c_double) + A = np.eye(2, dtype=c_double) + bupper = np.ones(2, dtype=c_double) + blower = -np.ones(2, dtype=c_double) + sense = np.zeros(2, dtype=c_int) + + for H in (H_symmetric, H_asymmetric, H_symmetric): + model = daqp.Model() + setup_flag, _ = model.setup( + H, f, A, bupper, blower, sense, is_avi=True) + self.assertEqual(setup_flag, 1) + x_model, _, flag_model, _ = model.solve() + x_ref, _, flag_ref, _ = daqp.solve( + H, f, A, bupper, blower, sense, is_avi=True) + self.assertEqual(flag_model, 1) + self.assertEqual(flag_ref, 1) + np.testing.assert_allclose(x_model, x_ref, atol=1.0e-8) + def test_warm_start_dual(self): """Dual warm start produces the same optimal solution as a cold start.""" H = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=c_double) diff --git a/src/api.c b/src/api.c index 70eadb0..05c30af 100644 --- a/src/api.c +++ b/src/api.c @@ -16,7 +16,7 @@ void daqp_solve(DAQPResult *res, DAQPWorkspace *work){ if(work->sing_ind != DAQP_UNCONSTRAINED_OPTIMAL){ // Select algorithm if(work->n_prox==0){ - if(work->avi == NULL){ + if(work->avi == NULL || work->avi->is_symmetric){ if(work->bnb != NULL) res->exitflag = daqp_bnb(work); else if(work->nh > 1) @@ -369,6 +369,8 @@ void allocate_daqp_ldp(DAQPWorkspace *work, int n, int m, int ms, int alloc_R, i #endif } void allocate_daqp_avi(DAQPAVI* avi, const int n){ + avi->is_symmetric = 0; + // Allocate matrices avi->Hsym = malloc(n*n*sizeof(c_float)); avi->Hs_rho = malloc(n*n*sizeof(c_float)); @@ -462,7 +464,9 @@ void daqp_extract_result(DAQPResult* res, DAQPWorkspace* work){ } // Shift back function value - if(work->v != NULL && work->avi == NULL && (work->Rinv != NULL || work->RinvD != NULL)){ // QP + if(work->v != NULL && + (work->avi == NULL || work->avi->is_symmetric) && + (work->Rinv != NULL || work->RinvD != NULL)){ // QP or symmetric AVI res->fval = work->fval; for(i=0;in;i++) res->fval-=work->v[i]*work->v[i]; res->fval *=0.5; diff --git a/src/eq_elim.c b/src/eq_elim.c index b0512e8..c8de3e6 100644 --- a/src/eq_elim.c +++ b/src/eq_elim.c @@ -125,7 +125,8 @@ static void apply_Q(const DAQPEqElim* eq, c_float* y){ * -------------------------------------------------------------------------*/ static int is_eq_elim_eligible(const DAQPWorkspace* work){ - if(work->avi != NULL || work->bnb != NULL || work->nh > 1) return 0; + if((work->avi != NULL && !work->avi->is_symmetric) || + work->bnb != NULL || work->nh > 1) return 0; if(work->qp == NULL || work->qp->A == NULL) return 0; /* * A singular Hessian is handled by the proximal method, which solves a diff --git a/src/utils.c b/src/utils.c index d6371ce..0cbb76d 100644 --- a/src/utils.c +++ b/src/utils.c @@ -58,12 +58,17 @@ int daqp_update_ldp(const int mask, DAQPWorkspace *work, DAQPProblem* qp){ if(work->avi == NULL) error_flag = daqp_update_Rinv(work, qp->H, qp->problem_type==2 ? 1 : 0); else{ - daqp_update_avi(work->avi,qp); - // Early unconstrained check for AVI: skip Cholesky if x=-H^{-1}f is feasible - int avi_unc = daqp_check_unconstrained(work,mask); - if(avi_unc == DAQP_UNCONSTRAINED_OPTIMAL) return 0; - daqp_lu(work->avi->H_rho, work->avi->P_H2, work->n); - error_flag = daqp_update_Rinv(work, work->avi->Hs_rho,0); + daqp_update_avi(work->avi,qp,work->settings->zero_tol); + if(work->avi->is_symmetric){ + error_flag = daqp_update_Rinv(work,qp->H,0); + } + else{ + // Early unconstrained check for AVI: skip Cholesky if x=-H^{-1}f is feasible + int avi_unc = daqp_check_unconstrained(work,mask); + if(avi_unc == DAQP_UNCONSTRAINED_OPTIMAL) return 0; + daqp_lu(work->avi->H_rho, work->avi->P_H2, work->n); + error_flag = daqp_update_Rinv(work, work->avi->Hs_rho,0); + } } if(error_flag<0) return error_flag; @@ -74,7 +79,8 @@ int daqp_update_ldp(const int mask, DAQPWorkspace *work, DAQPProblem* qp){ daqp_update_v(qp->f,work,mask); } - int unconstrained_flag = (work->avi != NULL) ? 1 : daqp_check_unconstrained(work,mask); + int unconstrained_flag = (work->avi != NULL && !work->avi->is_symmetric) + ? 1 : daqp_check_unconstrained(work,mask); if(unconstrained_flag == DAQP_UNCONSTRAINED_OPTIMAL) return 0; /* @@ -579,7 +585,7 @@ int daqp_check_unconstrained(DAQPWorkspace* work, const int mask){ // Compute x_unc stored temporarily in work->x. swp_ptr = work->x; work->u = work->xold; work->x = work->xold; work->xold = swp_ptr; - if(work->avi != NULL){ + if(work->avi != NULL && !work->avi->is_symmetric){ // AVI: unconstrained solution is x = -H^{-1} f if(work->qp->f != NULL) daqp_lu_solve(work->avi->LU_H, work->avi->P_H, work->qp->f, work->x, n); @@ -632,7 +638,7 @@ int daqp_check_unconstrained(DAQPWorkspace* work, const int mask){ return 1; } -int daqp_update_avi(DAQPAVI* avi, DAQPProblem* p){ +int daqp_update_avi(DAQPAVI* avi, DAQPProblem* p, c_float zero_tol){ const int n = p->n; // Setup matrices Hsym, Hs_rho, and H_rho, LU_H int i,j,disp; @@ -640,10 +646,15 @@ int daqp_update_avi(DAQPAVI* avi, DAQPProblem* p){ c_float min_diag = DAQP_INF; c_float max_row_sum = 0.0; c_float fro_norm_sq = 0.0; + c_float max_asymmetry = 0.0; avi->rho = 0.0; for (i = 0, disp=0; i < n; i++) { c_float row_sum = 0.0; for (j = 0; j < n; j++, disp++) { + if(j > i){ + c_float asymmetry = fabs(p->H[disp] - p->H[j * n + i]); + if(asymmetry > max_asymmetry) max_asymmetry = asymmetry; + } val = (p->H[disp] + p->H[j * n + i]) * 0.5; avi->Hsym[disp] = val; avi->Hs_rho[disp] = val; @@ -655,6 +666,11 @@ int daqp_update_avi(DAQPAVI* avi, DAQPProblem* p){ } if(row_sum > max_row_sum) max_row_sum = row_sum; } + c_float hessian_scale = sqrt(fro_norm_sq); + if(hessian_scale < 1.0) hessian_scale = 1.0; + avi->is_symmetric = max_asymmetry <= zero_tol * hessian_scale; + if(avi->is_symmetric) return 1; + // Regularization if(min_diag > 0.0 && max_row_sum > 0.0) avi->rho = sqrt(min_diag * max_row_sum);