Detect symmetric AVIs

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;

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);
 

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)

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;i<work->n;i++) res->fval-=work->v[i]*work->v[i];
         res->fval *=0.5;

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

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,18 +638,23 @@ 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;
     c_float val;
     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);