Line data    Source code

       1             : /**
       2             :  * @file linesearch.c
       3             :  * @author Ben Hermans
       4             :  * @brief Routines to perform exact or backtracking linesearch.
       5             :  * @details Once the direction is found using the semismooth Newton method, the functions in this file can 
       6             :  * be called to calculate a stepsize, either using exact linesearch or a backtracking linesearch to
       7             :  * satisfy the armijo condition.
       8             :  */
       9             : 
      10             : #include "linesearch.h"
      11             : #include "lin_alg.h"
      12             : #include <stdlib.h> //for sorting
      13             : 
      14         430 : c_float exact_linesearch(QPALMWorkspace *work, solver_common *c) {
      15         430 :     size_t n = work->data->n;
      16         430 :     size_t m = work->data->m;
      17             :     //Qd
      18         430 :     mat_vec(work->data->Q, work->solver->d, work->solver->Qd, c);
      19         430 :     if (work->settings->proximal) {
      20         322 :         vec_add_scaled(work->Qd, work->d, work->Qd, 1/work->gamma, n);
      21             :     }
      22             :     //Ad
      23         430 :     mat_vec(work->data->A, work->solver->d, work->solver->Ad, c);
      24             :     //eta = d'*Qd
      25         430 :     work->eta = vec_prod(work->d, work->Qd, n);
      26             :     //beta = d'*df
      27         430 :     work->beta = vec_prod(work->d, work->df, n);
      28             : 
      29             :     //delta = [-sqrt(sigma).*Ad; sqrt(sigma).*Ad]
      30         430 :     vec_ew_prod(work->sqrt_sigma, work->Ad, work->temp_m, m);
      31         430 :     prea_vec_copy(work->temp_m, work->delta + m, m); //shifted copy
      32         430 :     vec_self_mult_scalar(work->temp_m, -1, m);
      33         430 :     prea_vec_copy(work->temp_m, work->delta, m); 
      34             :     //alpha = [(y+sigma.*(Ax-bmin))./sigma_sqrt; (-y+sigma.*(bmax-Ax))./sigma_sqrt]
      35         430 :     vec_add_scaled(work->Ax, work->data->bmin, work->temp_m, -1, m);
      36         430 :     vec_ew_prod(work->sigma, work->temp_m, work->temp_m, m);
      37         430 :     vec_add_scaled(work->y, work->temp_m, work->temp_m, 1, m);
      38         430 :     vec_ew_div(work->temp_m, work->sqrt_sigma, work->temp_m, m);
      39         430 :     prea_vec_copy(work->temp_m, work->alpha, m);
      40         430 :     vec_add_scaled(work->data->bmax, work->Ax, work->temp_m, -1, m);
      41         430 :     vec_ew_prod(work->sigma, work->temp_m, work->temp_m, m);
      42         430 :     vec_add_scaled(work->temp_m, work->y, work->temp_m, -1, m);
      43         430 :     vec_ew_div(work->temp_m, work->sqrt_sigma, work->temp_m, m);
      44         430 :     prea_vec_copy(work->temp_m, work->alpha + m, m); //shifted copy
      45             :     // s = alpha./delta
      46         430 :     vec_ew_div(work->alpha, work->delta, work->temp_2m, m*2);
      47         430 :     vec_array_copy(work->temp_2m, work->s, m*2);
      48             :     // index_L = s>0
      49         430 :     size_t nL = 0;
      50             :     size_t i;
      51        4884 :     for (i = 0; i<m*2; i++){
      52        4454 :         if (work->temp_2m[i] > 0) {
      53        2883 :             work->index_L[i] = TRUE;
      54        2883 :             nL++;
      55             :         } else {
      56        1571 :             work->index_L[i] = FALSE;
      57             :         }     
      58             :     };
      59             : 
      60             :     //s = s(indL)
      61         430 :     select_subsequence(work->s, work->s, work->index_L, m*2);
      62             : 
      63             :     // index_P = delta > 0
      64        4884 :     for (i = 0; i<m*2; i++){
      65        4454 :         if (work->delta[i] > 0) {
      66        2017 :             work->index_P[i] = TRUE;
      67             :         } else {
      68        2437 :             work->index_P[i] = FALSE;
      69             :         }     
      70             :     };
      71             :     // index_J = (P&~L)|(~P&L);
      72        4884 :     for (i = 0; i<m*2; i++){
      73        4454 :         if ((work->index_P[i] + work->index_L[i]) == 1) {
      74        1142 :             work->index_J[i] = TRUE;
      75             :         } else {
      76        3312 :             work->index_J[i] = FALSE;
      77             :         }     
      78             :     };
      79             : 
      80             :     // a = eta+delta(J)'*delta(J);
      81             :     // b = beta-delta(J)'*alpha(J);
      82             :     c_float a, b;
      83         430 :     a = work->eta + vec_prod_ind(work->delta, work->delta, work->index_J, m*2);
      84         430 :     b = work->beta - vec_prod_ind(work->delta, work->alpha, work->index_J, m*2);
      85             :     // return 0;
      86             :     //s = sort(s)
      87             :     // qsort(work->s, m*2, sizeof(array_element), compare);
      88         430 :     qsort(work->s, nL, sizeof(array_element), compare);
      89             :     
      90         430 :     if (nL == 0 || a*work->s[0].x+b > 0) {
      91         243 :         return -b/a; 
      92             :     }
      93             : 
      94         187 :     i = 0;
      95             :     size_t iz;
      96         206 :     while (i < nL-1) {
      97         206 :         iz = work->s[i].i;
      98         206 :         if (work->index_P[iz]) {
      99         196 :             a = a + work->delta[iz]*work->delta[iz];
     100         196 :             b = b - work->delta[iz]*work->alpha[iz];
     101             :         } else {
     102          10 :             a = a - work->delta[iz]*work->delta[iz];
     103          10 :             b = b + work->delta[iz]*work->alpha[iz];
     104             :         }
     105         206 :         i++;
     106         206 :         if (a*work->s[i].x+b > 0) {
     107         187 :             return -b/a;
     108             :         }
     109             :     }
     110           0 :     iz = work->s[i].i;
     111           0 :     if (work->index_P[iz]) {
     112           0 :         a = a + work->delta[iz]*work->delta[iz];
     113           0 :         b = b - work->delta[iz]*work->alpha[iz];
     114             :     } else {
     115           0 :         a = a - work->delta[iz]*work->delta[iz];
     116           0 :         b = b + work->delta[iz]*work->alpha[iz];
     117             :     }
     118           0 :     return -b/a;
     119             : 
     120             : }
     121             : 
     122         430 : void vec_array_copy(c_float *a, array_element* b, size_t n) {
     123             :     size_t i;
     124             :     array_element ae;
     125             : 
     126        4884 :     for (i = 0; i < n; i++) {
     127        4454 :         ae.x = a[i];
     128        4454 :         ae.i = i;
     129        4454 :         b[i] = ae; 
     130             :     }
     131             : 
     132         430 : }
     133             : 
     134         430 : void select_subsequence(const array_element *a, array_element *b, const c_int *L, size_t n) {
     135             :     size_t i;
     136         430 :     size_t nb_elements = 0;
     137        4884 :     for (i = 0; i < n; i++) {
     138        4454 :         if (L[i]) {
     139        2883 :             b[nb_elements] = a[i];
     140        2883 :             nb_elements++;
     141             :         }
     142             :     }
     143         430 : }
     144             : 
     145         860 : c_float vec_prod_ind(const c_float *a, const c_float *b, const c_int *L, size_t n) {
     146         860 :   c_float prod = 0.0;
     147             :   size_t i;
     148             : 
     149        9768 :   for (i = 0; i < n; i++) {
     150        8908 :       if (L[i]) {
     151        2284 :           prod += a[i] * b[i];
     152             :       }   
     153             :   }
     154             : 
     155         860 :   return prod;
     156             : }
     157             : 
     158        5998 : int compare (const void * a, const void * b)
     159             : {
     160        5998 :     c_float f = ((struct array_element*)a)->x;
     161        5998 :     c_float s = ((struct array_element*)b)->x;
     162        5998 :     if (f > s) return  1;
     163        3248 :     if (f < s) return -1;
     164          38 :     return 0;
     165             : 
     166             : }