removed other norms and scaling_norm option from settings

laoyuer · Nov 6, 2017 · 1c2dae2 · 1c2dae2
1 parent 485abd0
commit 1c2dae2
Showing 21 changed files with 44 additions and 614 deletions.
diff --git a/docs/interfaces/solver_settings.rst b/docs/interfaces/solver_settings.rst
@@ -43,8 +43,6 @@ The solver settings are displayed in the following table. The settings marked wi
 +---------------------------------+---------------------------------------------------------------+--------------------------------------------------------------+-----------------------------------+
 | :code:`scaling`                 | Number of scaling iterations                                  | 0 (disabled) or 0 < :code:`scaling` (integer)                | 10                                |
 +---------------------------------+---------------------------------------------------------------+--------------------------------------------------------------+-----------------------------------+
-| :code:`scaling_norm`            | Scaling norm                                                  | 1, 2, -1 (infinity)                                          | -1 (infinity)                     |
-+---------------------------------+---------------------------------------------------------------+--------------------------------------------------------------+-----------------------------------+
 | :code:`adaptive_rho`            | Adaptive rho                                                  | True/False                                                   | True                              |
 +---------------------------------+---------------------------------------------------------------+--------------------------------------------------------------+-----------------------------------+
 | :code:`adaptive_rho_interval`   | Adaptive rho interval                                         | 0 (automatic) or 0 < :code:`adaptive_rho_interval` (integer) | 0                                 |

diff --git a/include/constants.h b/include/constants.h
@@ -70,8 +70,6 @@ static const char *LINSYS_SOLVER_NAME[] = {
 #define MAX_SCALING (1e+04)  ///< Maximum scaling value
 
 #if EMBEDDED != 1
-#define SCALING_NORM (-1)     ///< Scaling norm
-
 #define ADAPTIVE_RHO (1)
 #define ADAPTIVE_RHO_INTERVAL (0)
 #define ADAPTIVE_RHO_FRACTION (0.7)  ///< Fraction of setup time after which we update rho

diff --git a/include/lin_alg.h b/include/lin_alg.h
@@ -77,13 +77,6 @@ void vec_ew_max_vec(const c_float * a, const c_float * b, c_float * c, c_int n);
 /* Elementwise minimum between vectors c = min(a, b) */
 void vec_ew_min_vec(const c_float * a, const c_float * b, c_float * c, c_int n);
 
-/* Elementwise sqrt of sum of squares of vectors  */
-void vec_ew_sqrt_sos_vec(const c_float * a, const c_float * b,
-                         c_float * c, c_int n);
-
-/* Elementwise sum of vectors */
-void vec_ew_sum_vec(const c_float * a, const c_float * b,
-                    c_float * c, c_int n);
 #endif
 
 
@@ -102,30 +95,6 @@ i.e. scale the columns of A by d
 */
 void mat_postmult_diag(csc *A, const c_float *d);
 
-// #ifndef EMBEDDEED
-// TODO: Remove!
-/* Elementwise square matrix M */
-// void mat_ew_sq(csc * A);
-
-// TODO: Remove!
-/* Elementwise absolute value of matrix M */
-// void mat_ew_abs(csc * A);
-
-/**
- * Trace of matrix M in cdc format
- * @param  M Input matrix
- * @return   Trace
- */
-// c_float mat_trace(csc * M);
-
-/**
- * Frobenius norm squared of matrix M
- * @param  M Input matrix
- * @return   Frobenius norm squared
- */
-// c_float mat_fro_sq(csc * M);
-// #endif // ifndef embedded
-
 
 /* Matrix-vector multiplication
  *    y  =  A*x  (if plus_eq == 0)
@@ -173,57 +142,6 @@ void mat_inf_norm_rows(const csc * M, c_float * E);
 */
 void mat_inf_norm_cols_sym_triu(const csc * M, c_float * E);
 
-/**
-* 1 norm of each matrix column
-* @param M	Input matrix
-* @param E 	Vector of infinity norms
-*
-*/
-void mat_1_norm_cols(const csc * M, c_float * E);
-
-/**
-* 1 norm of each matrix row
-* @param M	Input matrix
-* @param E 	Vector of infinity norms
-*
-*/
-void mat_1_norm_rows(const csc * M, c_float * E);
-
-/**
-* 1 norm of each matrix column
-* Matrix M is symmetric upper-triangular
-*
-* @param M	Input matrix (symmetric, upper-triangular)
-* @param E 	Vector of infinity norms
-*
-*/
-void mat_1_norm_cols_sym_triu(const csc * M, c_float * E);
-
-/**
-* 2 norm of each matrix column
-* @param M	Input matrix
-* @param E 	Vector of infinity norms
-*
-*/
-void mat_2_norm_cols(const csc * M, c_float * E);
-
-/**
-* 2 norm of each matrix row
-* @param M	Input matrix
-* @param E 	Vector of infinity norms
-*
-*/
-void mat_2_norm_rows(const csc * M, c_float * E);
-
-/**
-* 2 norm of each matrix column
-* Matrix M is symmetric upper-triangular
-*
-* @param M	Input matrix (symmetric, upper-triangular)
-* @param E 	Vector of infinity norms
-*
-*/
-void mat_2_norm_cols_sym_triu(const csc * M, c_float * E);
 #endif  // EMBEDDED != 1
 
 /**

diff --git a/include/types.h b/include/types.h
@@ -139,7 +139,6 @@ typedef struct {
     c_int scaling;   ///< heuristic data scaling iterations. If 0, scaling disabled
 
 #if EMBEDDED != 1
-    c_int scaling_norm;  ///< scaling norm
     c_int adaptive_rho;  ///< boolean, is rho step size adaptive?
     c_int adaptive_rho_interval;  ///< Number of iterations between rho adaptations rho. If 0, it is automatic
     c_float adaptive_rho_tolerance;  ///< Tolerance X for adapting rho. The new rho has to be X times larger or 1/X times smaller than the current one to trigger a new factorization.

diff --git a/interfaces/matlab/codegen/render_workspace.m b/interfaces/matlab/codegen/render_workspace.m
@@ -71,7 +71,6 @@ function write_settings( f, settings, embedded_flag )
 fprintf(f, '%d, ',             settings.scaling);
 
 if embedded_flag ~= 1
-    fprintf(f, '%d, ', settings.scaling_norm);
     fprintf(f, '%d, ', settings.adaptive_rho);
     fprintf(f, '%d, ', settings.adaptive_rho_interval);
     fprintf(f, '(c_float)%.20f,', settings.adaptive_rho_tolerance);

diff --git a/interfaces/matlab/osqp.m b/interfaces/matlab/osqp.m
@@ -523,7 +523,7 @@ function codegen(this, target_dir, varargin)
 function currentSettings = validateSettings(this,isInitialization,varargin)
 
 %don't allow these fields to be changed
-unmodifiableFields = {'scaling', 'scaling_norm', 'linsys_solver'};
+unmodifiableFields = {'scaling', 'linsys_solver'};
 
 %get the current settings
 if(isInitialization)

diff --git a/interfaces/matlab/osqp_mex.cpp b/interfaces/matlab/osqp_mex.cpp
@@ -37,7 +37,6 @@ const char* OSQP_INFO_FIELDS[] = {"iter",         //c_int
 const char* OSQP_SETTINGS_FIELDS[] = {"rho",                        //c_float
                                       "sigma",                      //c_float
                                       "scaling",                    //c_int
-                                      "scaling_norm",               //c_int
                                       "adaptive_rho",               //c_int
                                       "adaptive_rho_interval",      //c_int
                                       "adaptive_rho_tolerance",     //c_float
@@ -789,7 +788,6 @@ mxArray* copySettingsToMxStruct(OSQPSettings* settings){
   mxSetField(mxPtr, 0, "rho",             mxCreateDoubleScalar(settings->rho));
   mxSetField(mxPtr, 0, "sigma",           mxCreateDoubleScalar(settings->sigma));
   mxSetField(mxPtr, 0, "scaling",         mxCreateDoubleScalar(settings->scaling));
-  mxSetField(mxPtr, 0, "scaling_norm",    mxCreateDoubleScalar(settings->scaling_norm));
   mxSetField(mxPtr, 0, "adaptive_rho",    mxCreateDoubleScalar(settings->adaptive_rho));
   mxSetField(mxPtr, 0, "adaptive_rho_interval",    mxCreateDoubleScalar(settings->adaptive_rho_interval));
   mxSetField(mxPtr, 0, "adaptive_rho_tolerance",    mxCreateDoubleScalar(settings->adaptive_rho_tolerance));
@@ -1054,7 +1052,6 @@ void copyMxStructToSettings(const mxArray* mxPtr, OSQPSettings* settings){
   settings->rho                       = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "rho"));
   settings->sigma                     = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "sigma"));
   settings->scaling                   = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "scaling"));
-  settings->scaling_norm              = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "scaling_norm"));
   settings->adaptive_rho              = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "adaptive_rho"));
   settings->adaptive_rho_interval     = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "adaptive_rho_interval"));
   settings->adaptive_rho_tolerance    = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "adaptive_rho_tolerance"));

diff --git a/interfaces/python/extension/include/osqpobjectpy.h b/interfaces/python/extension/include/osqpobjectpy.h
@@ -216,7 +216,7 @@ static PyObject * OSQP_setup(OSQP *self, PyObject *args, PyObject *kwargs) {
         static char *kwlist[] = {"dims",                          // nvars and ncons
                                  "Px", "Pi", "Pp", "q",           // Cost function
                                  "Ax", "Ai", "Ap", "l", "u",      // Constraints
-                                 "scaling", "scaling_norm",
+                                 "scaling", 
                                  "adaptive_rho", "adaptive_rho_interval",
 				 "adaptive_rho_tolerance", "adaptive_rho_fraction",
                                  "rho", "sigma", "max_iter", "eps_abs", "eps_rel", "eps_prim_inf", "eps_dual_inf", "alpha", "delta", "linsys_solver", "polish",
@@ -227,17 +227,17 @@ static PyObject * OSQP_setup(OSQP *self, PyObject *args, PyObject *kwargs) {
         #ifdef DLONG
 
         #ifdef DFLOAT
-        static char * argparse_string = "(LL)O!O!O!O!O!O!O!O!O!|LLLLffffLffffffLLLLLLL";
+        static char * argparse_string = "(LL)O!O!O!O!O!O!O!O!O!|LLLffffLffffffLLLLLLL";
         #else
-        static char * argparse_string = "(LL)O!O!O!O!O!O!O!O!O!|LLLLddddLddddddLLLLLLL";
+        static char * argparse_string = "(LL)O!O!O!O!O!O!O!O!O!|LLLddddLddddddLLLLLLL";
         #endif
 
         #else
 
         #ifdef DFLOAT
-        static char * argparse_string = "(ii)O!O!O!O!O!O!O!O!O!|iiiiffffiffffffiiiiiii";
+        static char * argparse_string = "(ii)O!O!O!O!O!O!O!O!O!|iiiffffiffffffiiiiiii";
         #else
-        static char * argparse_string = "(ii)O!O!O!O!O!O!O!O!O!|iiiiddddiddddddiiiiiii";
+        static char * argparse_string = "(ii)O!O!O!O!O!O!O!O!O!|iiiddddiddddddiiiiiii";
         #endif
 
         #endif
@@ -261,7 +261,6 @@ static PyObject * OSQP_setup(OSQP *self, PyObject *args, PyObject *kwargs) {
                                          &PyArray_Type, &l,
                                          &PyArray_Type, &u,
                                          &settings->scaling,
-                                         &settings->scaling_norm,
                                          &settings->adaptive_rho,
                                          &settings->adaptive_rho_interval,
                                          &settings->adaptive_rho_tolerance,

diff --git a/interfaces/python/extension/include/osqpworkspacepy.h b/interfaces/python/extension/include/osqpworkspacepy.h
@@ -214,11 +214,10 @@
      OSQPSettings *settings = self->workspace->settings;
 
      PyObject *return_dict = Py_BuildValue(
-         "{s:d,s:d,s:i,s:i,s:i,s:i,s:d,s:d,s:i,s:d,s:d,s:d, s:d, s:d, s:i, s:i, s:i, s:i}",
+         "{s:d,s:d,s:i,s:i,s:i,s:d,s:d,s:i,s:d,s:d,s:d, s:d, s:d, s:i, s:i, s:i, s:i}",
          "rho", (double)settings->rho,
          "sigma", (double)settings->sigma,
          "scaling", settings->scaling,
-         "scaling_norm", settings->scaling_norm,
          "adaptive_rho", settings->adaptive_rho,
          "adaptive_rho_interval", settings->adaptive_rho_interval,
          "adaptive_rho_tolerance", settings->adaptive_rho_tolerance,

diff --git a/interfaces/python/module/codegen/utils.py b/interfaces/python/module/codegen/utils.py
@@ -85,7 +85,6 @@ def write_settings(f, settings, name, embedded_flag):
     f.write("%d, " % settings['scaling'])
 
     if embedded_flag != 1:
-        f.write("%d, " % settings['scaling_norm'])
         f.write("%d, " % settings['adaptive_rho'])
         f.write("%d, " % settings['adaptive_rho_interval'])
         f.write("(c_float)%.20f, " % settings['adaptive_rho_tolerance'])

diff --git a/interfaces/python/module/interface.py b/interfaces/python/module/interface.py
@@ -124,25 +124,7 @@ def setup(self, P=None, q=None, A=None, l=None, u=None, **settings):
         l = np.maximum(l, -self._model.constant('OSQP_INFTY'))
 
         # Convert linsys_solver string to integer
-        linsys_solver_str = settings.pop('linsys_solver', '')
-        if not isinstance(linsys_solver_str, str):
-            raise TypeError("Setting linsys_solver " +
-                            "is required to be a string.")
-        linsys_solver_str = linsys_solver_str.lower()
-        if linsys_solver_str == 'suitesparse ldl':
-            settings['linsys_solver'] = \
-                    self._model.constant('SUITESPARSE_LDL_SOLVER')
-        elif linsys_solver_str == 'mkl pardiso':
-            settings['linsys_solver'] = self._model.constant('MKL_PARDISO_SOLVER')
-        # Default solver: Suitesparse LDL
-        elif linsys_solver_str == '':
-            settings['linsys_solver'] = \
-                    self._model.constant('SUITESPARSE_LDL_SOLVER')
-        else:   # default solver: Suitesparse LDL
-            warn("Linear system solver not recognized. " + 
-                 "Using default solver Suitesparse LDL.")
-            settings['linsys_solver'] = \
-                self._model.constant('SUITESPARSE_LDL_SOLVER')
+        self._linsys_solver_str_to_int(settings)
 
         self._model.setup((n, m), P.data, P.indices, P.indptr, q,
                           A.data, A.indices, A.indptr,
@@ -380,3 +362,25 @@ def codegen(self, folder, project_type='', parameters='vectors',
         # Generate code with codegen module
         cg.codegen(work, folder, python_ext_name, project_type,
                    embedded, force_rewrite, loop_unrolling)
+
+    def _linsys_solver_str_to_int(self, settings):
+        linsys_solver_str = settings.pop('linsys_solver', '')
+        if not isinstance(linsys_solver_str, str):
+            raise TypeError("Setting linsys_solver " +
+                            "is required to be a string.")
+        linsys_solver_str = linsys_solver_str.lower()
+        if linsys_solver_str == 'suitesparse ldl':
+            settings['linsys_solver'] = \
+                    self._model.constant('SUITESPARSE_LDL_SOLVER')
+        elif linsys_solver_str == 'mkl pardiso':
+            settings['linsys_solver'] = self._model.constant('MKL_PARDISO_SOLVER')
+        # Default solver: Suitesparse LDL
+        elif linsys_solver_str == '':
+            settings['linsys_solver'] = \
+                    self._model.constant('SUITESPARSE_LDL_SOLVER')
+        else:   # default solver: Suitesparse LDL
+            warn("Linear system solver not recognized. " + 
+                 "Using default solver Suitesparse LDL.")
+            settings['linsys_solver'] = \
+                self._model.constant('SUITESPARSE_LDL_SOLVER')
+
diff --git a/interfaces/python/modulepurepy/_osqp.py b/interfaces/python/modulepurepy/_osqp.py
@@ -127,7 +127,6 @@ class settings(object):
     -> These cannot be changed without running setup
     sigma    [1e-06]           - Regularization parameter for polish
     scaling  [10]            - Scaling/Equilibration iterations (0 disabled)
-    scaling_norm [-1]           - Equilibration scaling norm
 
     -> These can be changed without running setup
     rho  [1.6]                 - Step in ADMM procedure
@@ -151,7 +150,6 @@ def __init__(self, **kwargs):
         self.rho = kwargs.pop('rho', 0.1)
         self.sigma = kwargs.pop('sigma', 1e-06)
         self.scaling = kwargs.pop('scaling', 10)
-        self.scaling_norm = kwargs.pop('scaling_norm', -1)
         self.max_iter = kwargs.pop('max_iter', 4000)
         self.eps_abs = kwargs.pop('eps_abs', 1e-3)
         self.eps_rel = kwargs.pop('eps_rel', 1e-3)
@@ -337,34 +335,18 @@ def version(self):
         """
         return self._version
 
-    def _norm_KKT_cols(self, P, A, scaling_norm):
+    def _norm_KKT_cols(self, P, A):
         """
         Compute the norm of the KKT matrix from P and A
         """
 
-        if scaling_norm == 2:   # Scipy does not support sparse 2-norms
-            P = P.todense()
-            A = A.todense()
-            norm_function = np.linalg.norm
-        else:
-            norm_function = spspa.linalg.norm
-
-        if scaling_norm == -1:
-            scaling_norm = np.inf
-
         # First half
-        norm_P_cols = norm_function(P, scaling_norm, axis=0)
-        norm_A_cols = norm_function(A, scaling_norm, axis=0)
-        if scaling_norm == 2:
-            norm_first_half = np.sqrt(np.square(norm_P_cols) +
-                                      np.square(norm_A_cols))
-        elif scaling_norm == 1:
-            norm_first_half = norm_P_cols + norm_A_cols
-        elif scaling_norm == np.inf:  # Infinity norm
-            norm_first_half = np.maximum(norm_P_cols, norm_A_cols)
+        norm_P_cols = spspa.linalg.norm(P, np.inf, axis=0)
+        norm_A_cols = spspa.linalg.norm(A, np.inf, axis=0)
+        norm_first_half = np.maximum(norm_P_cols, norm_A_cols)
 
         # Second half (norm cols of A')
-        norm_second_half = norm_function(A, scaling_norm, axis=1)
+        norm_second_half = spspa.linalg.norm(A, np.inf, axis=1)
 
         return np.hstack((norm_first_half, norm_second_half))
 
@@ -401,17 +383,11 @@ def scale_data(self):
         """
         n = self.work.data.n
         m = self.work.data.m
-        scaling_norm = self.work.settings.scaling_norm
 
         # Initialize scaling
         s_temp = np.ones(n + m)
         c = 1.0  # Cost scaling
 
-        # Define reduced KKT matrix to scale
-        # KKT = spspa.vstack([
-        #       spspa.hstack([self.work.data.P, self.work.data.A.T]),
-        #       spspa.hstack([self.work.data.A,
-        #                     spspa.csc_matrix((m, m))])]).tocsc()
         # Define data
         P = self.work.data.P
         q = self.work.data.q
@@ -431,7 +407,7 @@ def scale_data(self):
         for i in range(self.work.settings.scaling):
 
             # First Step Ruiz
-            norm_cols = self._norm_KKT_cols(P, A, scaling_norm)
+            norm_cols = self._norm_KKT_cols(P, A)
             norm_cols = self._limit_scaling(norm_cols)  # Limit scaling
             sqrt_norm_cols = np.sqrt(norm_cols)         # Compute sqrt
             s_temp = np.reciprocal(sqrt_norm_cols)      # Elementwise recipr
@@ -602,11 +578,7 @@ def print_setup_header(self, data, settings):
               (settings.sigma, settings.alpha), end='')
         print("max_iter = %d" % settings.max_iter)
         if settings.scaling:
-            print("          scaling: on ", end='')
-            if settings.scaling_norm != -1:
-                print("(%d-norm), " % settings.scaling_norm, end='')
-            else:
-                print("(inf-norm), ", end='')
+            print("          scaling: on, ", end='')
         else:
             print("          scaling: off, ", end='')
         if settings.scaled_termination:

diff --git a/interfaces/python/tests/maros_meszaros/run_tests_osqp.py b/interfaces/python/tests/maros_meszaros/run_tests_osqp.py
@@ -213,7 +213,6 @@ def select_small_problems(problems):
             #  'scaled_termination': False,
             #  'adaptive_rho_interval': 100,
             #  'max_iter': 5000,
-            #  'scaling_norm': -1,
             'polish': False,
             #  'check_termination': 25
             }