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

SHOGUN Machine Learning Toolbox - Documentation