Bug fixes.

Author: wwkong

src/frameworks/IAIPAL.m

@@ -42,13 +42,12 @@
 
   % Initialize special constraint functions
   norm_fn = params.norm_fn;
-  pc_fn = params.set_projector; % projection onto the cone K.
 
   % Cone projector function on the point q = p + c * g(x) on the cones:
   %   -K (primal) and K^* (dual).
   function [dual_point, primal_point]= cone_proj(x, p, c)
     p_step = p + c * params.constr_fn(x);
-    primal_point = -pc_fn(-p_step);
+    primal_point = -params.set_projector(-p_step);
     dual_point = p_step - primal_point;
   end
 
@@ -65,7 +64,7 @@
       function val = wrap_f_s(x)
         [~, primal_proj_point]= cone_proj(x, p, c);
         p_step = p + c * params.constr_fn(x);
-        dist_val = norm_fn((-p_step) - primal_proj_point);
+        dist_val = norm_fn(p_step - primal_proj_point);
         val = 1 / (2 * c) * (dist_val ^ 2 - norm_fn(p) ^ 2);
       end
       function val = wrap_grad_f_s(x)

src/frameworks/penalty.m

@@ -191,6 +191,9 @@
   if (~isfield(params, 'i_logging')) 
     params.i_logging = false;
   end
+  if (~isfield(params, 'i_debug')) 
+    params.i_debug = false;
+  end
   if (~isfield(params, 'i_reset_prox_center')) 
     params.i_reset_prox_center = false;
   end

src/solvers/ACG.m

@@ -29,9 +29,8 @@
 %
 
   % Set some ACG global tolerances.
-  INEQ_COND_ERR_TOL = 1e-6;
+  INEQ_TOL = 1e-6;
   CURV_TOL = 1e-6;
-  DIV_TOL = 1e-6;
 
   %% PRE-PROCESSING
 
@@ -118,7 +117,7 @@
       error('Unknown ACG steptype!');
   end
 
-  % Safeguard against the case where (L == mu), i.e. nu = Inf.
+  % Safeguard against the case where (L == mu), i.e. lamK = Inf.
   L = max(L, mu + CURV_TOL);
 
   % Set up the oracle at x0
@@ -131,9 +130,9 @@
   function [local_L_est, aux_struct] = compute_approx_iter(L, mu, A_prev, y_prev, x_prev)
 
     % Simple quantities.
-    nu = nu_fn(mu, L);
-    nu_prev = (1 + mu * A_prev) * nu;
-    a_prev = (nu_prev + sqrt(nu_prev ^ 2 + 4 * nu_prev * A_prev)) / 2;
+    lamK = lamK_fn(mu, L);
+    tauK = (1 + mu * A_prev) * lamK;
+    a_prev = (tauK + sqrt(tauK ^ 2 + 4 * tauK * A_prev)) / 2;
     A = A_prev + a_prev;
     x_tilde_prev = (A_prev / A) * y_prev + a_prev / A * x_prev;
 
@@ -143,20 +142,26 @@
     grad_f_s_at_x_tilde_prev = o_x_tilde_prev.grad_f_s();
 
     % Oracle at y.
-    y_prox_mult = nu / (1 + nu * mu);
+    y_prox_mult = lamK / (1 + lamK * mu);
     y_prox_ctr = x_tilde_prev - y_prox_mult * grad_f_s_at_x_tilde_prev;
     [y, o_y] = get_y(y_prox_ctr, y_prox_mult);
     f_s_at_y = o_y.f_s();
 
     % Estimate of L based on y and x_tilde_prev.
-    local_L_est = max(0, 2 * (f_s_at_y - (f_s_at_x_tilde_prev + prod_fn(grad_f_s_at_x_tilde_prev, y - x_tilde_prev))) / ...
-                         norm_fn(y - x_tilde_prev) ^ 2);
+    LHS = f_s_at_y - (f_s_at_x_tilde_prev + prod_fn(grad_f_s_at_x_tilde_prev, y - x_tilde_prev));
+    dist_xt_y = norm_fn(y - x_tilde_prev);
+    RHS = L * dist_xt_y ^ 2 / 2;
+    local_L_est = max(0, 2 * LHS / dist_xt_y ^ 2);    
 
     % Save auxiliary quantities.
+    aux_struct.LHS = LHS;
+    aux_struct.RHS = RHS;
+    aux_struct.dist_xt_y = dist_xt_y;
+    aux_struct.descent_cond = (LHS <= RHS + INEQ_TOL);
     aux_struct.y = y;
     aux_struct.o_y = o_y;
-    aux_struct.nu = nu;
-    aux_struct.nu_prev = nu_prev;
+    aux_struct.lamK = lamK;
+    aux_struct.tauK = tauK;
     aux_struct.a_prev = a_prev;
     aux_struct.A = A;
     aux_struct.x_tilde_prev = x_tilde_prev;
@@ -206,22 +211,23 @@
       [local_L_est, aux_struct] = compute_approx_iter(L, mu, A_prev, y_prev, x_prev);
       iter = iter + 1;
 
-      % Adjust if the local estimate is smaller, up to a point.
-      if (L / 2 < local_L_est && local_L_est < L)
-        L = local_L_est;
-        [local_L_est, aux_struct] = compute_approx_iter(L, mu, A_prev, y_prev, x_prev);
-        iter = iter + 1;
-      end
-      
       % Update based on the value of the local L compared to the current estimate of L.
-      while (L < min([L_max, local_L_est]))
-        if (norm_fn(aux_struct.x_tilde_prev - aux_struct.y) ^ 2 <= DIV_TOL)
-          L = min(L_max, L * params.mult_L);
-        else
-          L = min(L_max, L * params.mult_L);
-        end
-        [local_L_est, aux_struct] = compute_approx_iter(L, mu, A_prev, y_prev, x_prev);
+      while (~aux_struct.descent_cond)        
+        L = min(L_max, L * params.mult_L);
+        [~, aux_struct] = compute_approx_iter(L, mu, A_prev, y_prev, x_prev);
         iter = iter + 1;
+        
+        % DEBUG ONLY
+        if (params.i_debug)
+          diff = abs(aux_struct.RHS - aux_struct.LHS);
+          dist_xt_y = aux_struct.dist_xt_y;
+          disp(table(local_L_est, L, L_max, dist_xt_y, diff, aux_struct.LHS, aux_struct.RHS , aux_struct.descent_cond));
+        end
+        % END DEBUG
+        
+        if (L >= L_max && ~aux_struct.descent_cond)
+          error('Theoretical upper bound on the upper curvature L_max does not appear to be correct!');
+        end
         if (toc(t_start) > time_limit)
           break;
         end
@@ -230,8 +236,8 @@
       % Load auxiliary quantities
       y = aux_struct.y;
       o_y = aux_struct.o_y;
-      nu = aux_struct.nu;
-      nu_prev = aux_struct.nu_prev;
+      lamK = aux_struct.lamK;
+      tauK = aux_struct.tauK;
       a_prev = aux_struct.a_prev;
       A = aux_struct.A;
       x_tilde_prev = aux_struct.x_tilde_prev;
@@ -242,9 +248,9 @@
     elseif strcmp(params.acg_steptype, "constant")
 
       % Iteration parameters.
-      nu = nu_fn(mu, L);
-      nu_prev = (1 + mu * A_prev) * nu;
-      a_prev = (nu_prev + sqrt(nu_prev ^ 2 + 4 * nu_prev * A_prev)) / 2;
+      lamK = lamK_fn(mu, L);
+      tauK = (1 + mu * A_prev) * lamK;
+      a_prev = (tauK + sqrt(tauK ^ 2 + 4 * tauK * A_prev)) / 2;
       A = A_prev + a_prev;
       x_tilde_prev = (A_prev / A) * y_prev + (a_prev / A) * x_prev;
 
@@ -254,7 +260,7 @@
       grad_f_s_at_x_tilde_prev = o_x_tilde_prev.grad_f_s();
 
       % Oracle at y.
-      y_prox_mult = nu / (1 + nu * mu);
+      y_prox_mult = lamK / (1 + lamK * mu);
       y_prox_ctr = x_tilde_prev - y_prox_mult * grad_f_s_at_x_tilde_prev;
       [y, o_y] = get_y(y_prox_ctr, y_prox_mult);
 
@@ -280,14 +286,14 @@
     %% COMPUTE (u, η), Γ, and x.
 
     % Compute x and u.
-    x = 1 / (1 + mu * A) * (x_prev  - a_prev / nu * (x_tilde_prev - y) + mu * (A_prev * x_prev + a_prev * y));
+    x = 1 / (1 + mu * A) * (x_prev  - a_prev / lamK * (x_tilde_prev - y) + mu * (A_prev * x_prev + a_prev * y));
     u = (x0 - x) / A;
 
     % Compute eta.
     if strcmp(params.eta_type, 'recursive')
       % Recursive
       gamma_at_x = f_n_at_y + f_s_at_x_tilde_prev + prod_fn(grad_f_s_at_x_tilde_prev, y - x_tilde_prev) + ...
-                   (mu / 2) * norm_fn(y - x_tilde_prev) ^ 2 + 1 / nu * prod_fn(x_tilde_prev - y, x - y) + ...
+                   (mu / 2) * norm_fn(y - x_tilde_prev) ^ 2 + 1 / lamK * prod_fn(x_tilde_prev - y, x - y) + ...
                    (mu / 2) * norm_fn(x - y) ^ 2;
       Gamma_at_x = a_prev / A * gamma_at_x + A_prev / A * Gamma_at_x_prev + A_prev / A * prod_fn(grad_Gamma_at_x_prev, x - x_prev) + ...
                    (mu / 2) * A_prev / A * norm_fn(x - x_prev) ^ 2;
@@ -296,8 +302,8 @@
     elseif strcmp(params.eta_type, 'accumulative')
       % Accumulative
       p_at_y = f_s_at_x_tilde_prev + prod_fn(grad_f_s_at_x_tilde_prev, y - x_tilde_prev) + mu / 2 * norm_fn(y - x_tilde_prev) ^ 2 + f_n_at_y;
-      sci = p_at_y + mu / 2 * norm_fn(y) ^ 2 - (1 / nu) * prod_fn(y, x_tilde_prev - y);
-      svi = - mu * y + (1 / nu) * (x_tilde_prev - y);
+      sci = p_at_y + mu / 2 * norm_fn(y) ^ 2 - (1 / lamK) * prod_fn(y, x_tilde_prev - y);
+      svi = - mu * y + (1 / lamK) * (x_tilde_prev - y);
       sni = mu / 2;
       scSum = scSum + a_prev * sci;
       svSum = svSum + a_prev * svi;
@@ -318,7 +324,7 @@
     % Check the negativity of eta in a relative sense.
     if strcmp(termination_type, "aipp")
       relative_exact_eta = exact_eta / max([norm_fn(u + x0 - y) ^ 2 / 2, 0.01]);
-      if (relative_exact_eta < -INEQ_COND_ERR_TOL)
+      if (relative_exact_eta < -INEQ_TOL)
         error(['eta is negative with a value of ', num2str(exact_eta)]);
       end
     end
@@ -334,7 +340,7 @@
       small_gd = f_at_y + prod_fn(u, x0 - y) - eta;
       del_gd = large_gd - small_gd;
       base = max([abs(large_gd), abs(small_gd), 0.01]);
-      if (del_gd / base < -INEQ_COND_ERR_TOL)
+      if (del_gd / base < -INEQ_TOL)
         model.status = -1;
         break;
       end
@@ -346,7 +352,7 @@
       large_gd = norm_fn(y - x0) ^ 2;
       del_gd = large_gd - small_gd;
       base = max([abs(large_gd), abs(small_gd), 0.01]);
-      if (del_gd / base < -INEQ_COND_ERR_TOL)
+      if (del_gd / base < -INEQ_TOL)
         model.status = -2;
         break;
       end
@@ -359,7 +365,7 @@
       large_gd = norm_fn(y - x0) ^ 2;
       del_gd = large_gd - small_gd;
       base = max([abs(large_gd), abs(small_gd), 0.01]);
-      if (del_gd / base < -INEQ_COND_ERR_TOL)
+      if (del_gd / base < -INEQ_TOL)
         model.status = -2;
         break;
       end
@@ -369,13 +375,13 @@
 
     % Termination for the AIPP method (Phase 1).
     if strcmp(termination_type, "aipp")
-      if (norm_fn(u) ^ 2 + 2 * eta <= sigma * norm_fn(x0 - y + u) ^ 2 + INEQ_COND_ERR_TOL)
+      if (norm_fn(u) ^ 2 + 2 * eta <= sigma * norm_fn(x0 - y + u) ^ 2 + INEQ_TOL)
         break;
       end
 
     % Termination for the AIPP method (with sigma square).
     elseif strcmp(termination_type, "aipp_sqr")
-      if (norm_fn(u) ^ 2 + 2 * eta <= sigma ^ 2 * norm_fn(x0 - y + u) ^ 2 + INEQ_COND_ERR_TOL)
+      if (norm_fn(u) ^ 2 + 2 * eta <= sigma ^ 2 * norm_fn(x0 - y + u) ^ 2 + INEQ_TOL)
         break;
       end
 
@@ -426,7 +432,7 @@
 
     % Update iterates.
     if (strcmp(params.eta_type, 'recursive'))
-      grad_gamma_at_x = 1 / nu * (x_tilde_prev - y) + mu * (x - y);
+      grad_gamma_at_x = 1 / lamK * (x_tilde_prev - y) + mu * (x - y);
       grad_Gamma_at_x = a_prev / A * grad_gamma_at_x + A_prev / A * grad_Gamma_at_x_prev + A_prev / A * mu * (x - x_prev);
       Gamma_at_x_prev = Gamma_at_x;
       grad_Gamma_at_x_prev = grad_Gamma_at_x;
@@ -466,8 +472,8 @@
 
 %% HELPER FUNCTIONS
 
-function out_nu = nu_fn(mu, L)
-  out_nu = 1 / (L - mu);
+function out_lamK =lamK_fn(mu, L)
+  out_lamK = 1 / (L - mu);
 end
 
 % Fills in parameters that were not set as input.

tests/papers/nl_iapial/extra/nl_iapial_qcqp_debug.m

@@ -0,0 +1,104 @@
+% Solve a multivariate nonconvex quadratically constrained quadratic programming  
+% problem constrained to a box using MULTIPLE SOLVERS.
+run('../../../../init.m');
+format long
+
+% Comment out later.
+dimN = 250;
+r = 1; 
+m = 1; 
+M = 1e4;
+
+% Run an instance via the command line.
+print_tbls(dimN, r, m, M);
+
+%% Utility functions
+function print_tbls(dimN, r, m, M) 
+
+  % Initialize
+  seed = 77777;
+  dimM = 10;
+  global_tol = 1e-5;
+  disp(table(dimN, r, m, M));
+  o_tbl = run_experiment(M, m, dimM, dimN, -r, r, seed, global_tol);
+  disp(o_tbl);
+  
+end
+function o_tbl = run_experiment(M, m, dimM, dimN, x_l, x_u, seed, global_tol)
+
+  [oracle, hparams] = test_fn_quad_box_constr_02(M, m, seed, dimM, dimN, x_l, x_u);
+  
+  % Set up the termination function.
+  function proj = proj_dh(a, b)
+    I1 = (abs(a - x_l) < 1e-12); 
+    I2 = (abs(a - x_u) <= 1e-12);
+    I3 = (abs(a - x_l) > 1e-12 & abs(a - x_u) > 1e-12);
+    proj = b;
+    proj(I1) = min(0, b(I1));
+    proj(I2) = max(0, b(I2));
+    proj(I3) = 0;
+  end
+  function proj = proj_NKt(~, b)
+    proj = min(0, b);
+  end
+  o_at_x0 = copy(oracle);
+  o_at_x0.eval(hparams.x0);
+  g0 = hparams.constr_fn(hparams.x0);
+  rho = global_tol * (1 + hparams.norm_fn(o_at_x0.grad_f_s()));
+  eta = global_tol * (1 + hparams.norm_fn(g0 - hparams.set_projector(g0)));
+  term_wrap = @(x,p) ...
+    termination_check(x, p, o_at_x0, hparams.constr_fn, hparams.grad_constr_fn, @proj_dh, @proj_NKt, hparams.norm_fn, rho, eta);
+
+  % Create the Model object and specify the solver.
+  ncvx_qc_qp = ConstrCompModel(oracle);
+  
+  % Set the curvatures and the starting point x0.
+  ncvx_qc_qp.x0 = hparams.x0;
+  ncvx_qc_qp.M = hparams.M;
+  ncvx_qc_qp.m = hparams.m;
+  ncvx_qc_qp.K_constr = hparams.K_constr;
+  ncvx_qc_qp.L_constr = hparams.L_constr;
+  
+  % Set the tolerances
+  ncvx_qc_qp.opt_tol = global_tol;
+  ncvx_qc_qp.feas_tol = global_tol;
+  ncvx_qc_qp.time_limit = 18000;
+  ncvx_qc_qp.iter_limit = 1000000;
+  
+  % Add linear constraints
+  ncvx_qc_qp.constr_fn = hparams.constr_fn;
+  ncvx_qc_qp.grad_constr_fn = hparams.grad_constr_fn;
+  ncvx_qc_qp.set_projector = hparams.set_projector;
+  ncvx_qc_qp.dual_cone_projector = hparams.dual_cone_projector;
+  
+  % Use a relative termination criterion.
+  ncvx_qc_qp.feas_type = 'relative';
+  ncvx_qc_qp.opt_type = 'relative';
+  
+  % Create some basic hparams.
+  base_hparam = struct();
+  base_hparam.i_debug = true;
+  base_hparam.termination_fn = term_wrap;
+  base_hparam.check_all_terminations = true;
+  
+  % Create the IAPIAL hparams.
+  ipl_hparam = base_hparam;
+  ipl_hparam.acg_steptype = 'constant';
+  ipl_hparam.sigma_min = 1/sqrt(2);
+  ipla_hparam = base_hparam;
+  ipla_hparam.acg_steptype = 'variable';
+  ipla_hparam.init_mult_L = 0.5;
+  ipla_hparam.sigma_min = 1/sqrt(2);
+  
+  % Run a benchmark test and print the summary.
+  hparam_arr = {ipla_hparam, ipl_hparam};
+  name_arr = {'IPL_A', 'IPL'};
+  framework_arr = {@IAIPAL, @IAIPAL};
+  solver_arr = {@ECG, @ECG};
+  
+  % Run the test.
+  [summary_tables, ~] = run_CCM_benchmark(ncvx_qc_qp, framework_arr, solver_arr, hparam_arr, name_arr);
+  o_tbl = summary_tables.all;
+
+end
+

tests/papers/nl_iapial/extra/nl_iapial_qcqp_iALM_only_params.m

@@ -56,7 +56,7 @@ function print_tbls(dimN, r, m, M)
   % Set the tolerances
   ncvx_qc_qp.opt_tol = global_tol;
   ncvx_qc_qp.feas_tol = global_tol;
-  ncvx_qc_qp.time_limit = Inf;
+  ncvx_qc_qp.time_limit = 6000;
   ncvx_qc_qp.iter_limit = 1000000;
 
   % Add linear constraints