SHOGUN  6.0.0
altra.cpp
Go to the documentation of this file.
1 /* This program is free software: you can redistribute it and/or modify
3  * the Free Software Foundation, either version 3 of the License, or
4  * (at your option) any later version.
5  *
6  * This program is distributed in the hope that it will be useful,
7  * but WITHOUT ANY WARRANTY; without even the implied warranty of
8  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
9  * GNU General Public License for more details.
10  *
11  * You should have received a copy of the GNU General Public License
12  * along with this program. If not, see <http://www.gnu.org/licenses/>.
13  *
14  * Copyright (C) 2009 - 2012 Jun Liu and Jieping Ye
15  */
16
17
19 #ifdef USE_GPL_SHOGUN
21 #include <stdlib.h>
22 #include <string.h>
23
24 void altra(double *x, double *v, int n, double *ind, int nodes, double mult)
25 {
26  int i, j;
27  double lambda,twoNorm, ratio;
28
29  /*
30  * test whether the first node is special
31  */
32  if ((int) ind[0]==-1){
33
34  /*
35  *Recheck whether ind[1] equals to zero
36  */
37  if ((int) ind[1]!=-1){
38  printf("\n Error! \n Check ind");
39  exit(1);
40  }
41
42  lambda=mult*ind[2];
43
44  for(j=0;j<n;j++){
45  if (v[j]>lambda)
46  x[j]=v[j]-lambda;
47  else
48  if (v[j]<-lambda)
49  x[j]=v[j]+lambda;
50  else
51  x[j]=0;
52  }
53
54  i=1;
55  }
56  else{
57  shogun::sg_memcpy(x, v, sizeof(double) * n);
58  i=0;
59  }
60
61  /*
62  * sequentially process each node
63  *
64  */
65  for(;i < nodes; i++){
66  /*
67  * compute the L2 norm of this group
68  */
69  twoNorm=0;
70  for(j=(int) ind[3*i]-1;j< (int) ind[3*i+1];j++)
71  twoNorm += x[j] * x[j];
72  twoNorm=sqrt(twoNorm);
73
74  lambda=mult*ind[3*i+2];
75  if (twoNorm>lambda){
76  ratio=(twoNorm-lambda)/twoNorm;
77
78  /*
79  * shrinkage this group by ratio
80  */
81  for(j=(int) ind[3*i]-1;j<(int) ind[3*i+1];j++)
82  x[j]*=ratio;
83  }
84  else{
85  /*
86  * threshold this group to zero
87  */
88  for(j=(int) ind[3*i]-1;j<(int) ind[3*i+1];j++)
89  x[j]=0;
90  }
91  }
92 }
93
94 void altra_mt(double *X, double *V, int n, int k, double *ind, int nodes, double mult)
95 {
96  int i, j;
97
98  double *x=(double *)malloc(sizeof(double)*k);
99  double *v=(double *)malloc(sizeof(double)*k);
100
101  for (i=0;i<n;i++){
102  /*
103  * copy a row of V to v
104  *
105  */
106  for(j=0;j<k;j++)
107  v[j]=V[j*n + i];
108
109  altra(x, v, k, ind, nodes, mult);
110
111  /*
112  * copy the solution to X
113  */
114  for(j=0;j<k;j++)
115  X[j*n+i]=x[j];
116  }
117
118  free(x);
119  free(v);
120 }
121
122 void computeLambda2Max(double *lambda2_max, double *x, int n, double *ind, int nodes)
123 {
124  int i, j;
125  double twoNorm;
126
127  *lambda2_max=0;
128
129  for(i=0;i < nodes; i++){
130  /*
131  * compute the L2 norm of this group
132  */
133  twoNorm=0;
134  for(j=(int) ind[3*i]-1;j< (int) ind[3*i+1];j++)
135  twoNorm += x[j] * x[j];
136  twoNorm=sqrt(twoNorm);
137
138  twoNorm=twoNorm/ind[3*i+2];
139
140  if (twoNorm >*lambda2_max )
141  *lambda2_max=twoNorm;
142  }
143 }
144
145 double treeNorm(double *x, int ldx, int n, double *ind, int nodes){
146
147  int i, j;
148  double twoNorm, lambda;
149
150  double tree_norm = 0;
151
152  /*
153  * test whether the first node is special
154  */
155  if ((int) ind[0]==-1){
156
157  /*
158  *Recheck whether ind[1] equals to zero
159  */
160  if ((int) ind[1]!=-1){
161  printf("\n Error! \n Check ind");
162  exit(1);
163  }
164
165  lambda=ind[2];
166
167  for(j=0;j<n*ldx;j+=ldx){
168  tree_norm+=fabs(x[j]);
169  }
170
171  tree_norm = tree_norm * lambda;
172
173  i=1;
174  }
175  else{
176  i=0;
177  }
178
179  /*
180  * sequentially process each node
181  *
182  */
183  for(;i < nodes; i++){
184  /*
185  * compute the L2 norm of this group
186  */
187  twoNorm=0;
188
189  int n_in_node = (int) ind[3*i+1] - (int) ind[3*i]-1;
190  for(j=(int) ind[3*i]-1;j< (int) ind[3*i]-1 + n_in_node*ldx;j+=ldx)
191  twoNorm += x[j] * x[j];
192  twoNorm=sqrt(twoNorm);
193
194  lambda=ind[3*i+2];
195
196  tree_norm = tree_norm + lambda*twoNorm;
197  }
198
199  return tree_norm;
200 }
201
202 double findLambdaMax(double *v, int n, double *ind, int nodes){
203
204  int i;
205  double lambda=0,squaredWeight=0, lambda1,lambda2;
206  double *x=(double *)malloc(sizeof(double)*n);
207  double *ind2=(double *)malloc(sizeof(double)*nodes*3);
208  int num=0;
209
210  for(i=0;i<n;i++){
211  lambda+=v[i]*v[i];
212  }
213
214  if ( (int)ind[0]==-1 )
215  squaredWeight=n*ind[2]*ind[2];
216  else
217  squaredWeight=ind[2]*ind[2];
218
219  for (i=1;i<nodes;i++){
220  squaredWeight+=ind[3*i+2]*ind[3*i+2];
221  }
222
223  /* set lambda to an initial guess
224  */
225  lambda=sqrt(lambda/squaredWeight);
226
227  /*
228  printf("\n\n lambda=%2.5f",lambda);
229  */
230
231  /*
232  *copy ind to ind2,
233  *and scale the weight 3*i+2
234  */
235  for(i=0;i<nodes;i++){
236  ind2[3*i]=ind[3*i];
237  ind2[3*i+1]=ind[3*i+1];
238  ind2[3*i+2]=ind[3*i+2]*lambda;
239  }
240
241  /* test whether the solution is zero or not
242  */
243  altra(x, v, n, ind2, nodes);
244  for(i=0;i<n;i++){
245  if (x[i]!=0)
246  break;
247  }
248
249  if (i>=n) {
250  /*x is a zero vector*/
251  lambda2=lambda;
252  lambda1=lambda;
253
254  num=0;
255
256  while(1){
257  num++;
258
259  lambda2=lambda;
260  lambda1=lambda1/2;
261  /* update ind2
262  */
263  for(i=0;i<nodes;i++){
264  ind2[3*i+2]=ind[3*i+2]*lambda1;
265  }
266
267  /* compute and test whether x is zero
268  */
269  altra(x, v, n, ind2, nodes);
270  for(i=0;i<n;i++){
271  if (x[i]!=0)
272  break;
273  }
274
275  if (i<n){
276  break;
277  /*x is not zero
278  *we have found lambda1
279  */
280  }
281  }
282
283  }
284  else{
285  /*x is a non-zero vector*/
286  lambda2=lambda;
287  lambda1=lambda;
288
289  num=0;
290  while(1){
291  num++;
292
293  lambda1=lambda2;
294  lambda2=lambda2*2;
295  /* update ind2
296  */
297  for(i=0;i<nodes;i++){
298  ind2[3*i+2]=ind[3*i+2]*lambda2;
299  }
300
301  /* compute and test whether x is zero
302  */
303  altra(x, v, n, ind2, nodes);
304  for(i=0;i<n;i++){
305  if (x[i]!=0)
306  break;
307  }
308
309  if (i>=n){
310  break;
311  /*x is a zero vector
312  *we have found lambda2
313  */
314  }
315  }
316  }
317
318  /*
319  printf("\n num=%d, lambda1=%2.5f, lambda2=%2.5f",num, lambda1,lambda2);
320  */
321
322  while ( fabs(lambda2-lambda1) > lambda2 * 1e-10 ){
323
324  num++;
325
326  lambda=(lambda1+lambda2)/2;
327
328  /* update ind2
329  */
330  for(i=0;i<nodes;i++){
331  ind2[3*i+2]=ind[3*i+2]*lambda;
332  }
333
334  /* compute and test whether x is zero
335  */
336  altra(x, v, n, ind2, nodes);
337  for(i=0;i<n;i++){
338  if (x[i]!=0)
339  break;
340  }
341
342  if (i>=n){
343  lambda2=lambda;
344  }
345  else{
346  lambda1=lambda;
347  }
348
349  /*
350  printf("\n lambda1=%2.5f, lambda2=%2.5f",lambda1,lambda2);
351  */
352  }
353
354  /*
355  printf("\n num=%d",num);
356
357  printf(" lambda1=%2.5f, lambda2=%2.5f",lambda1,lambda2);
358
359 */
360
361  free(x);
362  free(ind2);
363
364  return lambda2;
365 }
366
367 double findLambdaMax_mt(double *V, int n, int k, double *ind, int nodes)
368 {
369  int i, j;
370
371  double *v=(double *)malloc(sizeof(double)*k);
372  double lambda;
373
374  double lambdaMax=0;
375
376  for (i=0;i<n;i++){
377  /*
378  * copy a row of V to v
379  *
380  */
381  for(j=0;j<k;j++)
382  v[j]=V[j*n + i];
383
384  lambda = findLambdaMax(v, k, ind, nodes);
385
386  /*
387  printf("\n lambda=%5.2f",lambda);
388  */
389
390  if (lambda>lambdaMax)
391  lambdaMax=lambda;
392  }
393
394  /*
395  printf("\n *lambdaMax=%5.2f",*lambdaMax);
396  */
397
398  free(v);
399  return lambdaMax;
400 }
401
402 #endif //USE_GPL_SHOGUN

SHOGUN Machine Learning Toolbox - Documentation