source: pro-violet-viettel/sourcecode/application/libraries/PHPExcel/Shared/JAMA/QRDecomposition.php @ 769

Last change on this file since 769 was 289, checked in by dungnv, 11 years ago
File size: 5.7 KB
RevLine 
[289]1<?php
2/**
3 *      @package JAMA
4 *
5 *      For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
6 *      orthogonal matrix Q and an n-by-n upper triangular matrix R so that
7 *      A = Q*R.
8 *
9 *      The QR decompostion always exists, even if the matrix does not have
10 *      full rank, so the constructor will never fail.  The primary use of the
11 *      QR decomposition is in the least squares solution of nonsquare systems
12 *      of simultaneous linear equations.  This will fail if isFullRank()
13 *      returns false.
14 *
15 *      @author  Paul Meagher
16 *      @license PHP v3.0
17 *      @version 1.1
18 */
19class PHPExcel_Shared_JAMA_QRDecomposition {
20
21        const MatrixRankException       = "Can only perform operation on full-rank matrix.";
22
23        /**
24         *      Array for internal storage of decomposition.
25         *      @var array
26         */
27        private $QR = array();
28
29        /**
30         *      Row dimension.
31         *      @var integer
32         */
33        private $m;
34
35        /**
36        *       Column dimension.
37        *       @var integer
38        */
39        private $n;
40
41        /**
42         *      Array for internal storage of diagonal of R.
43         *      @var  array
44         */
45        private $Rdiag = array();
46
47
48        /**
49         *      QR Decomposition computed by Householder reflections.
50         *
51         *      @param matrix $A Rectangular matrix
52         *      @return Structure to access R and the Householder vectors and compute Q.
53         */
54        public function __construct($A) {
55                if($A instanceof PHPExcel_Shared_JAMA_Matrix) {
56                        // Initialize.
57                        $this->QR = $A->getArrayCopy();
58                        $this->m  = $A->getRowDimension();
59                        $this->n  = $A->getColumnDimension();
60                        // Main loop.
61                        for ($k = 0; $k < $this->n; ++$k) {
62                                // Compute 2-norm of k-th column without under/overflow.
63                                $nrm = 0.0;
64                                for ($i = $k; $i < $this->m; ++$i) {
65                                        $nrm = hypo($nrm, $this->QR[$i][$k]);
66                                }
67                                if ($nrm != 0.0) {
68                                        // Form k-th Householder vector.
69                                        if ($this->QR[$k][$k] < 0) {
70                                                $nrm = -$nrm;
71                                        }
72                                        for ($i = $k; $i < $this->m; ++$i) {
73                                                $this->QR[$i][$k] /= $nrm;
74                                        }
75                                        $this->QR[$k][$k] += 1.0;
76                                        // Apply transformation to remaining columns.
77                                        for ($j = $k+1; $j < $this->n; ++$j) {
78                                                $s = 0.0;
79                                                for ($i = $k; $i < $this->m; ++$i) {
80                                                        $s += $this->QR[$i][$k] * $this->QR[$i][$j];
81                                                }
82                                                $s = -$s/$this->QR[$k][$k];
83                                                for ($i = $k; $i < $this->m; ++$i) {
84                                                        $this->QR[$i][$j] += $s * $this->QR[$i][$k];
85                                                }
86                                        }
87                                }
88                                $this->Rdiag[$k] = -$nrm;
89                        }
90                } else {
91                        throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
92                }
93        }       //      function __construct()
94
95
96        /**
97         *      Is the matrix full rank?
98         *
99         *      @return boolean true if R, and hence A, has full rank, else false.
100         */
101        public function isFullRank() {
102                for ($j = 0; $j < $this->n; ++$j) {
103                        if ($this->Rdiag[$j] == 0) {
104                                return false;
105                        }
106                }
107                return true;
108        }       //      function isFullRank()
109
110
111        /**
112         *      Return the Householder vectors
113         *
114         *      @return Matrix Lower trapezoidal matrix whose columns define the reflections
115         */
116        public function getH() {
117                for ($i = 0; $i < $this->m; ++$i) {
118                        for ($j = 0; $j < $this->n; ++$j) {
119                                if ($i >= $j) {
120                                        $H[$i][$j] = $this->QR[$i][$j];
121                                } else {
122                                        $H[$i][$j] = 0.0;
123                                }
124                        }
125                }
126                return new PHPExcel_Shared_JAMA_Matrix($H);
127        }       //      function getH()
128
129
130        /**
131         *      Return the upper triangular factor
132         *
133         *      @return Matrix upper triangular factor
134         */
135        public function getR() {
136                for ($i = 0; $i < $this->n; ++$i) {
137                        for ($j = 0; $j < $this->n; ++$j) {
138                                if ($i < $j) {
139                                        $R[$i][$j] = $this->QR[$i][$j];
140                                } elseif ($i == $j) {
141                                        $R[$i][$j] = $this->Rdiag[$i];
142                                } else {
143                                        $R[$i][$j] = 0.0;
144                                }
145                        }
146                }
147                return new PHPExcel_Shared_JAMA_Matrix($R);
148        }       //      function getR()
149
150
151        /**
152         *      Generate and return the (economy-sized) orthogonal factor
153         *
154         *      @return Matrix orthogonal factor
155         */
156        public function getQ() {
157                for ($k = $this->n-1; $k >= 0; --$k) {
158                        for ($i = 0; $i < $this->m; ++$i) {
159                                $Q[$i][$k] = 0.0;
160                        }
161                        $Q[$k][$k] = 1.0;
162                        for ($j = $k; $j < $this->n; ++$j) {
163                                if ($this->QR[$k][$k] != 0) {
164                                        $s = 0.0;
165                                        for ($i = $k; $i < $this->m; ++$i) {
166                                                $s += $this->QR[$i][$k] * $Q[$i][$j];
167                                        }
168                                        $s = -$s/$this->QR[$k][$k];
169                                        for ($i = $k; $i < $this->m; ++$i) {
170                                                $Q[$i][$j] += $s * $this->QR[$i][$k];
171                                        }
172                                }
173                        }
174                }
175                /*
176                for($i = 0; $i < count($Q); ++$i) {
177                        for($j = 0; $j < count($Q); ++$j) {
178                                if(! isset($Q[$i][$j]) ) {
179                                        $Q[$i][$j] = 0;
180                                }
181                        }
182                }
183                */
184                return new PHPExcel_Shared_JAMA_Matrix($Q);
185        }       //      function getQ()
186
187
188        /**
189         *      Least squares solution of A*X = B
190         *
191         *      @param Matrix $B A Matrix with as many rows as A and any number of columns.
192         *      @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
193         */
194        public function solve($B) {
195                if ($B->getRowDimension() == $this->m) {
196                        if ($this->isFullRank()) {
197                                // Copy right hand side
198                                $nx = $B->getColumnDimension();
199                                $X  = $B->getArrayCopy();
200                                // Compute Y = transpose(Q)*B
201                                for ($k = 0; $k < $this->n; ++$k) {
202                                        for ($j = 0; $j < $nx; ++$j) {
203                                                $s = 0.0;
204                                                for ($i = $k; $i < $this->m; ++$i) {
205                                                        $s += $this->QR[$i][$k] * $X[$i][$j];
206                                                }
207                                                $s = -$s/$this->QR[$k][$k];
208                                                for ($i = $k; $i < $this->m; ++$i) {
209                                                        $X[$i][$j] += $s * $this->QR[$i][$k];
210                                                }
211                                        }
212                                }
213                                // Solve R*X = Y;
214                                for ($k = $this->n-1; $k >= 0; --$k) {
215                                        for ($j = 0; $j < $nx; ++$j) {
216                                                $X[$k][$j] /= $this->Rdiag[$k];
217                                        }
218                                        for ($i = 0; $i < $k; ++$i) {
219                                                for ($j = 0; $j < $nx; ++$j) {
220                                                        $X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
221                                                }
222                                        }
223                                }
224                                $X = new PHPExcel_Shared_JAMA_Matrix($X);
225                                return ($X->getMatrix(0, $this->n-1, 0, $nx));
226                        } else {
227                                throw new PHPExcel_Calculation_Exception(self::MatrixRankException);
228                        }
229                } else {
230                        throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
231                }
232        }       //      function solve()
233
234}       //      PHPExcel_Shared_JAMA_class QRDecomposition
Note: See TracBrowser for help on using the repository browser.