← Back to branch summary

~ubuntu-branches/ubuntu/lucid/phpmyadmin/lucid

« back to all changes in this revision

Viewing changes to libraries/PHPExcel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php

  • Committer: Bazaar Package Importer
  • Author(s): Michal Čihař
  • Date: 2010-03-08 15:25:00 UTC
  • mfrom: (1.2.8 upstream)
  • Revision ID: james.westby@ubuntu.com-20100308152500-6e8hmuqc5co39de5
Tags: 4:3.3.0-1
* New upstream version.
* Rediff debian/patches.
* Fix permissions on mediawiki export extension.

Show diffs side-by-side

added added

removed removed

Lines of Context:
libraries/PHPExcel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php
 
1
<?php
 
2
/**
 
3
* @package JAMA
 
4
*
 
5
* For an m-by-n matrix A with m >= n, the singular value decomposition is
 
6
* an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
 
7
* an n-by-n orthogonal matrix V so that A = U*S*V'.
 
8
*
 
9
* The singular values, sigma[$k] = S[$k][$k], are ordered so that
 
10
* sigma[0] >= sigma[1] >= ... >= sigma[n-1].
 
11
*
 
12
* The singular value decompostion always exists, so the constructor will
 
13
* never fail.  The matrix condition number and the effective numerical
 
14
* rank can be computed from this decomposition.
 
15
*
 
16
* @author  Paul Meagher
 
17
* @license PHP v3.0
 
18
* @version 1.1
 
19
*/
 
20
class SingularValueDecomposition  {
 
21
 
 
22
  /**
 
23
  * Internal storage of U.
 
24
  * @var array
 
25
  */
 
26
  var $U = array();
 
27
 
 
28
  /**
 
29
  * Internal storage of V.
 
30
  * @var array
 
31
  */
 
32
  var $V = array();
 
33
 
 
34
  /**
 
35
  * Internal storage of singular values.
 
36
  * @var array
 
37
  */
 
38
  var $s = array();
 
39
 
 
40
  /**
 
41
  * Row dimension.
 
42
  * @var int
 
43
  */
 
44
  var $m;
 
45
 
 
46
  /**
 
47
  * Column dimension.
 
48
  * @var int
 
49
  */
 
50
  var $n;
 
51
 
 
52
  /**
 
53
  * Construct the singular value decomposition
 
54
  *
 
55
  * Derived from LINPACK code.
 
56
  *
 
57
  * @param $A Rectangular matrix
 
58
  * @return Structure to access U, S and V.
 
59
  */
 
60
  function SingularValueDecomposition ($Arg) {
 
61
 
 
62
    // Initialize.
 
63
 
 
64
    $A = $Arg->getArrayCopy();
 
65
    $this->m = $Arg->getRowDimension();
 
66
    $this->n = $Arg->getColumnDimension();
 
67
    $nu      = min($this->m, $this->n);
 
68
    $e       = array();
 
69
    $work    = array();
 
70
    $wantu   = true;
 
71
    $wantv   = true;      
 
72
    $nct = min($this->m - 1, $this->n);
 
73
    $nrt = max(0, min($this->n - 2, $this->m));   
 
74
 
 
75
    // Reduce A to bidiagonal form, storing the diagonal elements
 
76
    // in s and the super-diagonal elements in e.
 
77
 
 
78
    for ($k = 0; $k < max($nct,$nrt); $k++) {
 
79
 
 
80
      if ($k < $nct) {
 
81
        // Compute the transformation for the k-th column and
 
82
        // place the k-th diagonal in s[$k].
 
83
        // Compute 2-norm of k-th column without under/overflow.
 
84
        $this->s[$k] = 0;
 
85
        for ($i = $k; $i < $this->m; $i++)
 
86
          $this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
 
87
        if ($this->s[$k] != 0.0) {
 
88
          if ($A[$k][$k] < 0.0)
 
89
            $this->s[$k] = -$this->s[$k];
 
90
          for ($i = $k; $i < $this->m; $i++)
 
91
            $A[$i][$k] /= $this->s[$k];
 
92
          $A[$k][$k] += 1.0;
 
93
        }
 
94
        $this->s[$k] = -$this->s[$k];
 
95
      }           
 
96
 
 
97
      for ($j = $k + 1; $j < $this->n; $j++) {
 
98
        if (($k < $nct) & ($this->s[$k] != 0.0)) {
 
99
          // Apply the transformation.
 
100
          $t = 0;
 
101
          for ($i = $k; $i < $this->m; $i++)
 
102
            $t += $A[$i][$k] * $A[$i][$j];
 
103
          $t = -$t / $A[$k][$k];
 
104
          for ($i = $k; $i < $this->m; $i++)
 
105
            $A[$i][$j] += $t * $A[$i][$k];
 
106
          // Place the k-th row of A into e for the
 
107
          // subsequent calculation of the row transformation.
 
108
          $e[$j] = $A[$k][$j];
 
109
        }
 
110
      }
 
111
 
 
112
      if ($wantu AND ($k < $nct)) {
 
113
        // Place the transformation in U for subsequent back
 
114
        // multiplication.
 
115
        for ($i = $k; $i < $this->m; $i++)
 
116
          $this->U[$i][$k] = $A[$i][$k];
 
117
      }
 
118
 
 
119
      if ($k < $nrt) {
 
120
        // Compute the k-th row transformation and place the
 
121
        // k-th super-diagonal in e[$k].
 
122
        // Compute 2-norm without under/overflow.
 
123
        $e[$k] = 0;
 
124
        for ($i = $k + 1; $i < $this->n; $i++)
 
125
          $e[$k] = hypo($e[$k], $e[$i]);
 
126
        if ($e[$k] != 0.0) {
 
127
          if ($e[$k+1] < 0.0)
 
128
            $e[$k] = -$e[$k];
 
129
          for ($i = $k + 1; $i < $this->n; $i++)
 
130
            $e[$i] /= $e[$k];
 
131
          $e[$k+1] += 1.0;
 
132
        }
 
133
        $e[$k] = -$e[$k];
 
134
        if (($k+1 < $this->m) AND ($e[$k] != 0.0)) {
 
135
          // Apply the transformation.
 
136
          for ($i = $k+1; $i < $this->m; $i++)
 
137
            $work[$i] = 0.0;
 
138
          for ($j = $k+1; $j < $this->n; $j++)
 
139
            for ($i = $k+1; $i < $this->m; $i++)
 
140
              $work[$i] += $e[$j] * $A[$i][$j];
 
141
          for ($j = $k + 1; $j < $this->n; $j++) {
 
142
            $t = -$e[$j] / $e[$k+1];
 
143
            for ($i = $k + 1; $i < $this->m; $i++)
 
144
              $A[$i][$j] += $t * $work[$i];
 
145
          }
 
146
        }
 
147
        if ($wantv) {
 
148
          // Place the transformation in V for subsequent
 
149
          // back multiplication.
 
150
          for ($i = $k + 1; $i < $this->n; $i++)
 
151
            $this->V[$i][$k] = $e[$i];
 
152
        }        
 
153
      }
 
154
    } 
 
155
       
 
156
    // Set up the final bidiagonal matrix or order p.
 
157
    $p = min($this->n, $this->m + 1);
 
158
    if ($nct < $this->n)
 
159
      $this->s[$nct] = $A[$nct][$nct];
 
160
    if ($this->m < $p)
 
161
      $this->s[$p-1] = 0.0;
 
162
    if ($nrt + 1 < $p)
 
163
      $e[$nrt] = $A[$nrt][$p-1];
 
164
    $e[$p-1] = 0.0;
 
165
    // If required, generate U.
 
166
    if ($wantu) {
 
167
      for ($j = $nct; $j < $nu; $j++) {
 
168
        for ($i = 0; $i < $this->m; $i++)
 
169
          $this->U[$i][$j] = 0.0;
 
170
        $this->U[$j][$j] = 1.0;
 
171
      }
 
172
      for ($k = $nct - 1; $k >= 0; $k--) {
 
173
        if ($this->s[$k] != 0.0) {
 
174
          for ($j = $k + 1; $j < $nu; $j++) {
 
175
            $t = 0;
 
176
            for ($i = $k; $i < $this->m; $i++)
 
177
              $t += $this->U[$i][$k] * $this->U[$i][$j];
 
178
            $t = -$t / $this->U[$k][$k];
 
179
            for ($i = $k; $i < $this->m; $i++)
 
180
              $this->U[$i][$j] += $t * $this->U[$i][$k];
 
181
          }
 
182
          for ($i = $k; $i < $this->m; $i++ )
 
183
            $this->U[$i][$k] = -$this->U[$i][$k];
 
184
          $this->U[$k][$k] = 1.0 + $this->U[$k][$k];
 
185
          for ($i = 0; $i < $k - 1; $i++)
 
186
             $this->U[$i][$k] = 0.0;
 
187
         } else {
 
188
           for ($i = 0; $i < $this->m; $i++)
 
189
             $this->U[$i][$k] = 0.0;
 
190
           $this->U[$k][$k] = 1.0;
 
191
         }
 
192
      }
 
193
    }
 
194
 
 
195
    // If required, generate V.
 
196
    if ($wantv) {
 
197
      for ($k = $this->n - 1; $k >= 0; $k--) {
 
198
        if (($k < $nrt) AND ($e[$k] != 0.0)) {
 
199
          for ($j = $k + 1; $j < $nu; $j++) {
 
200
            $t = 0;
 
201
            for ($i = $k + 1; $i < $this->n; $i++)
 
202
              $t += $this->V[$i][$k]* $this->V[$i][$j];
 
203
            $t = -$t / $this->V[$k+1][$k];
 
204
            for ($i = $k + 1; $i < $this->n; $i++)
 
205
              $this->V[$i][$j] += $t * $this->V[$i][$k];
 
206
          }
 
207
        }
 
208
        for ($i = 0; $i < $this->n; $i++)
 
209
           $this->V[$i][$k] = 0.0;
 
210
        $this->V[$k][$k] = 1.0;
 
211
      }
 
212
    }
 
213
        
 
214
    // Main iteration loop for the singular values.
 
215
    $pp   = $p - 1;
 
216
    $iter = 0;
 
217
    $eps  = pow(2.0, -52.0);
 
218
    while ($p > 0) {      
 
219
      
 
220
      // Here is where a test for too many iterations would go.
 
221
      // This section of the program inspects for negligible
 
222
      // elements in the s and e arrays.  On completion the
 
223
      // variables kase and k are set as follows:
 
224
      // kase = 1  if s(p) and e[k-1] are negligible and k<p
 
225
      // kase = 2  if s(k) is negligible and k<p
 
226
      // kase = 3  if e[k-1] is negligible, k<p, and
 
227
      //           s(k), ..., s(p) are not negligible (qr step).
 
228
      // kase = 4  if e(p-1) is negligible (convergence).
 
229
 
 
230
      for ($k = $p - 2; $k >= -1; $k--) {
 
231
        if ($k == -1)
 
232
          break;
 
233
        if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) {
 
234
          $e[$k] = 0.0;
 
235
          break;
 
236
        }
 
237
      }
 
238
      if ($k == $p - 2)
 
239
        $kase = 4;
 
240
      else {
 
241
        for ($ks = $p - 1; $ks >= $k; $ks--) {
 
242
          if ($ks == $k)
 
243
            break;
 
244
          $t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.);
 
245
          if (abs($this->s[$ks]) <= $eps * $t)  {
 
246
            $this->s[$ks] = 0.0;
 
247
            break;
 
248
          }
 
249
        }
 
250
        if ($ks == $k)
 
251
           $kase = 3;
 
252
        else if ($ks == $p-1)
 
253
           $kase = 1;
 
254
        else {
 
255
           $kase = 2;
 
256
           $k = $ks;
 
257
        }
 
258
      }
 
259
      $k++;
 
260
 
 
261
      // Perform the task indicated by kase.
 
262
      switch ($kase) {
 
263
         // Deflate negligible s(p).
 
264
         case 1:
 
265
           $f = $e[$p-2];
 
266
           $e[$p-2] = 0.0;
 
267
           for ($j = $p - 2; $j >= $k; $j--) {
 
268
             $t  = hypo($this->s[$j],$f);
 
269
             $cs = $this->s[$j] / $t;
 
270
             $sn = $f / $t;
 
271
             $this->s[$j] = $t;
 
272
             if ($j != $k) {
 
273
               $f = -$sn * $e[$j-1];
 
274
               $e[$j-1] = $cs * $e[$j-1];
 
275
             }
 
276
             if ($wantv) {
 
277
               for ($i = 0; $i < $this->n; $i++) {
 
278
                 $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1];
 
279
                 $this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1];
 
280
                 $this->V[$i][$j] = $t;
 
281
               }
 
282
             }
 
283
           }
 
284
           break;
 
285
         // Split at negligible s(k).
 
286
         case 2:
 
287
           $f = $e[$k-1];
 
288
           $e[$k-1] = 0.0;
 
289
           for ($j = $k; $j < $p; $j++) {
 
290
             $t = hypo($this->s[$j], $f);
 
291
             $cs = $this->s[$j] / $t;
 
292
             $sn = $f / $t;
 
293
             $this->s[$j] = $t;
 
294
             $f = -$sn * $e[$j];
 
295
             $e[$j] = $cs * $e[$j];
 
296
             if ($wantu) {
 
297
               for ($i = 0; $i < $this->m; $i++) {
 
298
                 $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1];
 
299
                 $this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1];
 
300
                 $this->U[$i][$j] = $t;
 
301
               }
 
302
             }
 
303
           }
 
304
           break;           
 
305
         // Perform one qr step.
 
306
         case 3:                  
 
307
           // Calculate the shift.                          
 
308
           $scale = max(max(max(max(
 
309
                    abs($this->s[$p-1]),abs($this->s[$p-2])),abs($e[$p-2])),
 
310
                    abs($this->s[$k])), abs($e[$k]));
 
311
           $sp   = $this->s[$p-1] / $scale;
 
312
           $spm1 = $this->s[$p-2] / $scale;
 
313
           $epm1 = $e[$p-2] / $scale;
 
314
           $sk   = $this->s[$k] / $scale;
 
315
           $ek   = $e[$k] / $scale;
 
316
           $b    = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
 
317
           $c    = ($sp * $epm1) * ($sp * $epm1);
 
318
           $shift = 0.0;
 
319
           if (($b != 0.0) || ($c != 0.0)) {
 
320
             $shift = sqrt($b * $b + $c);
 
321
             if ($b < 0.0)
 
322
               $shift = -$shift;
 
323
             $shift = $c / ($b + $shift);
 
324
           }
 
325
           $f = ($sk + $sp) * ($sk - $sp) + $shift;
 
326
           $g = $sk * $ek;  
 
327
           // Chase zeros.  
 
328
           for ($j = $k; $j < $p-1; $j++) {
 
329
             $t  = hypo($f,$g);
 
330
             $cs = $f/$t;
 
331
             $sn = $g/$t;
 
332
             if ($j != $k)
 
333
               $e[$j-1] = $t;
 
334
             $f = $cs * $this->s[$j] + $sn * $e[$j];
 
335
             $e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
 
336
             $g = $sn * $this->s[$j+1];
 
337
             $this->s[$j+1] = $cs * $this->s[$j+1];
 
338
             if ($wantv) {
 
339
               for ($i = 0; $i < $this->n; $i++) {
 
340
                 $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1];
 
341
                 $this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1];
 
342
                 $this->V[$i][$j] = $t;
 
343
               }
 
344
             }                                                                      
 
345
             $t = hypo($f,$g);
 
346
             $cs = $f/$t;
 
347
             $sn = $g/$t;
 
348
             $this->s[$j] = $t;
 
349
             $f = $cs * $e[$j] + $sn * $this->s[$j+1];
 
350
             $this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1];
 
351
             $g = $sn * $e[$j+1];
 
352
             $e[$j+1] = $cs * $e[$j+1];
 
353
             if ($wantu && ($j < $this->m - 1)) {
 
354
               for ($i = 0; $i < $this->m; $i++) {
 
355
                 $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1];
 
356
                 $this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1];
 
357
                 $this->U[$i][$j] = $t;
 
358
               }
 
359
             }                
 
360
           }
 
361
           $e[$p-2] = $f;
 
362
           $iter = $iter + 1;
 
363
           break;
 
364
         // Convergence.
 
365
         case 4:
 
366
           // Make the singular values positive.
 
367
           if ($this->s[$k] <= 0.0) {
 
368
             $this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
 
369
             if ($wantv) {
 
370
               for ($i = 0; $i <= $pp; $i++)
 
371
                 $this->V[$i][$k] = -$this->V[$i][$k];
 
372
             }
 
373
           }
 
374
           // Order the singular values.  
 
375
           while ($k < $pp) {
 
376
            if ($this->s[$k] >= $this->s[$k+1])
 
377
              break;
 
378
            $t = $this->s[$k];
 
379
            $this->s[$k] = $this->s[$k+1];
 
380
            $this->s[$k+1] = $t;
 
381
            if ($wantv AND ($k < $this->n - 1)) {
 
382
              for ($i = 0; $i < $this->n; $i++) {
 
383
                $t = $this->V[$i][$k+1];
 
384
                $this->V[$i][$k+1] = $this->V[$i][$k];
 
385
                $this->V[$i][$k] = $t;
 
386
              }
 
387
            }
 
388
            if ($wantu AND ($k < $this->m-1)) {
 
389
              for ($i = 0; $i < $this->m; $i++) {
 
390
                $t = $this->U[$i][$k+1];
 
391
                $this->U[$i][$k+1] = $this->U[$i][$k];
 
392
                $this->U[$i][$k] = $t;
 
393
              }
 
394
            }
 
395
            $k++;
 
396
           }
 
397
           $iter = 0;
 
398
           $p--;
 
399
           break;                                   
 
400
      } // end switch      
 
401
    } // end while
 
402
 
 
403
    /*
 
404
    echo "<p>Output A</p>";   
 
405
    $A = new Matrix($A);        
 
406
    $A->toHTML();
 
407
    
 
408
    echo "<p>Matrix U</p>";    
 
409
    echo "<pre>";
 
410
    print_r($this->U);        
 
411
    echo "</pre>";
 
412
    
 
413
    echo "<p>Matrix V</p>";    
 
414
    echo "<pre>";
 
415
    print_r($this->V);        
 
416
    echo "</pre>";    
 
417
    
 
418
    echo "<p>Vector S</p>";    
 
419
    echo "<pre>";
 
420
    print_r($this->s);        
 
421
    echo "</pre>";    
 
422
    exit;
 
423
    */
 
424
    
 
425
  } // end constructor
 
426
  
 
427
  /**
 
428
  * Return the left singular vectors
 
429
  * @access public
 
430
  * @return U
 
431
  */
 
432
  function getU() {
 
433
    return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
 
434
  }
 
435
 
 
436
  /**
 
437
  * Return the right singular vectors
 
438
  * @access public   
 
439
  * @return V
 
440
  */
 
441
  function getV() {
 
442
    return new Matrix($this->V);
 
443
  }
 
444
 
 
445
  /**
 
446
  * Return the one-dimensional array of singular values
 
447
  * @access public   
 
448
  * @return diagonal of S.
 
449
  */
 
450
  function getSingularValues() {
 
451
    return $this->s;
 
452
  }
 
453
 
 
454
  /**
 
455
  * Return the diagonal matrix of singular values
 
456
  * @access public   
 
457
  * @return S
 
458
  */
 
459
  function getS() {
 
460
    for ($i = 0; $i < $this->n; $i++) {
 
461
      for ($j = 0; $j < $this->n; $j++)
 
462
        $S[$i][$j] = 0.0;
 
463
      $S[$i][$i] = $this->s[$i];
 
464
    }
 
465
    return new Matrix($S);
 
466
  }
 
467
 
 
468
  /**
 
469
  * Two norm
 
470
  * @access public   
 
471
  * @return max(S)
 
472
  */
 
473
  function norm2() {
 
474
    return $this->s[0];
 
475
  }
 
476
 
 
477
  /**
 
478
  * Two norm condition number
 
479
  * @access public   
 
480
  * @return max(S)/min(S)
 
481
  */
 
482
  function cond() {
 
483
    return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
 
484
  }
 
485
 
 
486
  /**
 
487
  * Effective numerical matrix rank
 
488
  * @access public   
 
489
  * @return Number of nonnegligible singular values.
 
490
  */
 
491
  function rank() {
 
492
    $eps = pow(2.0, -52.0);
 
493
    $tol = max($this->m, $this->n) * $this->s[0] * $eps;
 
494
    $r = 0;
 
495
    for ($i = 0; $i < count($this->s); $i++) {
 
496
      if ($this->s[$i] > $tol)
 
497
        $r++;
 
498
    }
 
499
    return $r;
 
500
  }
 
501
}