GestureRecognitionToolkit  Version: 1.0 Revision: 04-03-15
The Gesture Recognition Toolkit (GRT) is a cross-platform, open-source, c++ machine learning library for real-time gesture recognition.
SVD.cpp
1 /*
2  GRT MIT License
3  Copyright (c) <2012> <Nicholas Gillian, Media Lab, MIT>
4 
5  Permission is hereby granted, free of charge, to any person obtaining a copy of this software
6  and associated documentation files (the "Software"), to deal in the Software without restriction,
7  including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense,
8  and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so,
9  subject to the following conditions:
10 
11  The above copyright notice and this permission notice shall be included in all copies or substantial
12  portions of the Software.
13 
14  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT
15  LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
16  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
17  WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
18  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
19  */
20 #include "SVD.h"
21 
22 namespace GRT{
23 
24 bool SVD::solve(Matrix<double> &a){
25 
26  //Setup the memory
27  m = a.getNumRows();
28  n = a.getNumCols();
29  u = a;
30  v.resize(n,n);
31  w.resize(n);
32 
33  eps = numeric_limits< double >::epsilon();
34  if( !decompose() ) return false;
35  if( !reorder() ) return false;
36 
37  tsh = 0.5*sqrt(m+n+1.)*w[0]*eps;
38 
39  return true;
40 }
41 
42 bool SVD::solveVector(vector <double> &b, vector <double> &x, double thresh) {
43  UINT i,j,jj;
44  double s;
45  if(b.size() != m || x.size() != n){
46  return false;
47  }
48  vector <double> tmp(n);
49  tsh = (thresh >= 0. ? thresh : 0.5*sqrt(m+n+1.)*w[0]*eps);
50  for (j=0;j<n;j++) {
51  s=0.0;
52  if (w[j] > tsh) {
53  for (i=0;i<m;i++) s += u[i][j]*b[i];
54  s /= w[j];
55  }
56  tmp[j]=s;
57  }
58  for (j=0;j<n;j++) {
59  s=0.0;
60  for (jj=0;jj<n;jj++) s += v[j][jj]*tmp[jj];
61  x[j]=s;
62  }
63  return true;
64 }
65 
66 bool SVD::solve(Matrix <double> &b, Matrix <double> &x, double thresh){
67  UINT i,j,m=b.getNumCols();
68  if (b.getNumRows() != n || x.getNumRows() != n || b.getNumCols() != x.getNumCols()){
69  return false;
70  }
71  vector <double> xx(n);
72  for (j=0;j<m;j++) {
73  for (i=0;i<n;i++) xx[i] = b[i][j];
74  solveVector(xx,xx,thresh);
75  for (i=0;i<n;i++) x[i][j] = xx[i];
76  }
77  return true;
78 }
79 UINT SVD::rank(double thresh) {
80  UINT j,nr=0;
81  tsh = (thresh >= 0. ? thresh : 0.5*sqrt(m+n+1.)*w[0]*eps);
82  for (j=0;j<n;j++) if (w[j] > tsh) nr++;
83  return nr;
84 }
85 
86 UINT SVD::nullity(double thresh) {
87  UINT j,nn=0;
88  tsh = (thresh >= 0. ? thresh : 0.5*sqrt(m+n+1.)*w[0]*eps);
89  for (j=0;j<n;j++) if (w[j] <= tsh) nn++;
90  return nn;
91 }
92 
93 Matrix <double> SVD::range(double thresh){
94  UINT i,j,nr=0;
95  Matrix <double> rnge(m,rank(thresh));
96  for (j=0;j<n;j++) {
97  if (w[j] > tsh) {
98  for (i=0;i<m;i++) rnge[i][nr] = u[i][j];
99  nr++;
100  }
101  }
102  return rnge;
103 }
104 
105 Matrix <double> SVD::nullspace(double thresh){
106  UINT j,jj,nn=0;
107  Matrix <double> nullsp(n,nullity(thresh));
108  for (j=0;j<n;j++) {
109  if (w[j] <= tsh) {
110  for (jj=0;jj<n;jj++) nullsp[jj][nn] = v[jj][j];
111  nn++;
112  }
113  }
114  return nullsp;
115 }
116 
117 double SVD::inv_condition() {
118  return (w[0] <= 0. || w[n-1] <= 0.) ? 0. : w[n-1]/w[0];
119 }
120 
121 bool SVD::decompose() {
122  bool flag;
123  int i,its,j,jj,k,l,nm,N,M;
124  double anorm,c,f,g,h,s,scale,x,y,z;
125  vector <double> rv1(n);
126  g = scale = anorm = 0.0;
127  N = int(n);
128  M = int(m);
129  l = 0;
130  for (i=0;i<N;i++) {
131  l=i+2;
132  rv1[i]=scale*g;
133  g=s=scale=0.0;
134  if (i < M) {
135  for (k=i;k<M;k++) scale += fabs(u[k][i]);
136  if (scale != 0.0) {
137  for (k=i;k<M;k++) {
138  u[k][i] /= scale;
139  s+= u[k][i]*u[k][i];
140  }
141  f=u[i][i];
142  g = -SIGN(sqrt(s),f);
143  h=f*g-s;
144  u[i][i]=f-g;
145  for (j=l-1;j<N;j++) {
146  for (s=0.0,k=i;k<M;k++) s += u[k][i] * u[k][j];
147  f=s/h;
148  for (k=i;k<M;k++) u[k][j] += f * u[k][i];
149  }
150  for (k=i;k<M;k++) u[k][i] *= scale;
151  }
152  }
153  w[i]=scale *g;
154  g=s=scale=0.0;
155  if (i+1 <= M && i+1 != N) {
156  for (k=l-1;k<N;k++) scale += fabs(u[i][k]);
157  if (scale != 0.0) {
158  for (k=l-1;k<N;k++) {
159  u[i][k] /= scale;
160  s += u[i][k]*u[i][k];
161  }
162  f=u[i][l-1];
163  g = -SIGN(sqrt(s),f);
164  h=f*g-s;
165  u[i][l-1]=f-g;
166  for (k=l-1;k<N;k++) rv1[k]=u[i][k]/h;
167  for (j=l-1;j<M;j++) {
168  for (s=0.0,k=l-1;k<N;k++) s += u[j][k]*u[i][k];
169  for (k=l-1;k<N;k++) u[j][k] += s*rv1[k];
170  }
171  for (k=l-1;k<N;k++) u[i][k]*= scale;
172  }
173  }
174  anorm=MAX(anorm,(fabs(w[i])+fabs(rv1[i])));
175  }
176  for (i=N-1;i>=0;i--) {
177  if (i < N-1) {
178  if (g != 0.0) {
179  for (j=l;j<N;j++)
180  v[j][i]=u[i][j]/u[i][l]/g;
181  for (j=l;j<N;j++) {
182  for (s=0.0,k=l;k<N;k++) s += u[i][k]*v[k][j];
183  for (k=l;k<N;k++) v[k][j] += s*v[k][i];
184  }
185  }
186  for (j=l;j<N;j++) v[i][j]=v[j][i]=0.0;
187  }
188  v[i][i]=1.0;
189  g=rv1[i];
190  l=i;
191  }
192  for (i=MIN(M,N)-1;i>=0;i--) {
193  l=i+1;
194  g=w[i];
195  for (j=l;j<N;j++) u[i][j]=0.0;
196  if (g != 0.0) {
197  g=1.0/g;
198  for (j=l;j<N;j++) {
199  for (s=0.0,k=l;k<M;k++) s += u[k][i]*u[k][j];
200  f=(s/u[i][i])*g;
201  for (k=i;k<M;k++) u[k][j] += f*u[k][i];
202  }
203  for (j=i;j<M;j++) u[j][i] *= g;
204  } else for (j=i;j<M;j++) u[j][i] =0.0;
205  ++u[i][i];
206  }
207  for (k=N-1;k>=0;k--) {
208  for (its=0;its<MAX_NUM_SVD_ITER;its++) {
209  flag=true;
210  for (l=k;l>=0;l--) {
211  nm=l-1;
212  if (l == 0 || fabs(rv1[l]) <= eps*anorm) {
213  flag=false;
214  break;
215  }
216  if (fabs(w[nm]) <= eps*anorm) break;
217  }
218  if (flag) {
219  c=0.0;
220  s=1.0;
221  for (i=l;i<k+1;i++) {
222  f=s*rv1[i];
223  rv1[i]=c*rv1[i];
224  if (fabs(f) <= eps*anorm) break;
225  g=w[i];
226  h=pythag(f,g);
227  w[i]=h;
228  h=1.0/h;
229  c=g*h;
230  s = -f*h;
231  for (j=0;j<M;j++) {
232  y=u[j][nm];
233  z=u[j][i];
234  u[j][nm]=y*c+z*s;
235  u[j][i]=z*c-y*s;
236  }
237  }
238  }
239  z=w[k];
240  if (l == k) {
241  if (z < 0.0) {
242  w[k] = -z;
243  for (j=0;j<N;j++) v[j][k] = -v[j][k];
244  }
245  break;
246  }
247  if (its == MAX_NUM_SVD_ITER-1){
248  return false;
249  }
250  x=w[l];
251  nm=k-1;
252  y=w[nm];
253  g=rv1[nm];
254  h=rv1[k];
255  f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
256  g=pythag(f,1.0);
257  f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
258  c=s=1.0;
259  for (j=l;j<=nm;j++) {
260  i=j+1;
261  g=rv1[i];
262  y=w[i];
263  h=s*g;
264  g=c*g;
265  z=pythag(f,h);
266  rv1[j]=z;
267  c=f/z;
268  s=h/z;
269  f=x*c+g*s;
270  g=g*c-x*s;
271  h=y*s;
272  y *= c;
273  for (jj=0;jj<N;jj++) {
274  x=v[jj][j];
275  z=v[jj][i];
276  v[jj][j]=x*c+z*s;
277  v[jj][i]=z*c-x*s;
278  }
279  z=pythag(f,h);
280  w[j]=z;
281  if (z) {
282  z=1.0/z;
283  c=f*z;
284  s=h*z;
285  }
286  f=c*g+s*y;
287  x=c*y-s*g;
288  for (jj=0;jj<M;jj++) {
289  y=u[jj][j];
290  z=u[jj][i];
291  u[jj][j]=y*c+z*s;
292  u[jj][i]=z*c-y*s;
293  }
294  }
295  rv1[l]=0.0;
296  rv1[k]=f;
297  w[k]=x;
298  }
299  }
300 
301  return true;
302 }
303 
304 bool SVD::reorder() {
305  UINT i,j,k,s,inc=1;
306  double sw;
307  vector <double> su(m);
308  vector <double> sv(n);
309  do { inc *= 3; inc++; } while (inc <= n);
310  do {
311  inc /= 3;
312  for (i=inc;i<n;i++) {
313  sw = w[i];
314  for (k=0;k<m;k++) su[k] = u[k][i];
315  for (k=0;k<n;k++) sv[k] = v[k][i];
316  j = i;
317  while (w[j-inc] < sw) {
318  w[j] = w[j-inc];
319  for (k=0;k<m;k++) u[k][j] = u[k][j-inc];
320  for (k=0;k<n;k++) v[k][j] = v[k][j-inc];
321  j -= inc;
322  if (j < inc) break;
323  }
324  w[j] = sw;
325  for (k=0;k<m;k++) u[k][j] = su[k];
326  for (k=0;k<n;k++) v[k][j] = sv[k];
327  }
328  } while (inc > 1);
329  for (k=0;k<n;k++) {
330  s=0;
331  for (i=0;i<m;i++) if (u[i][k] < 0.) s++;
332  for (j=0;j<n;j++) if (v[j][k] < 0.) s++;
333  if (s > (m+n)/2) {
334  for (i=0;i<m;i++) u[i][k] = -u[i][k];
335  for (j=0;j<n;j++) v[j][k] = -v[j][k];
336  }
337  }
338  return true;
339 }
340 
341 double SVD::pythag(const double a, const double b) {
342  double absa=fabs(a);
343  double absb=fabs(b);
344  return (absa > absb ? absa*sqrt(1.0+SQR(absb/absa)) : (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb))));
345 }
346 
347 }
Definition: AdaBoost.cpp:25
virtual bool resize(const unsigned int r, const unsigned int c)
Definition: Matrix.h:234