Proyectos de Subversion Moodle

Rev

| Ultima modificación | Ver Log |

Rev Autor Línea Nro. Línea
1 efrain 1 
<?php
2 
 
3 
namespace Matrix\Decomposition;
4 
 
5 
use Matrix\Exception;
6 
use Matrix\Matrix;
7 
 
8 
class QR
9 
{
10 
    private $qrMatrix;
11 
    private $rows;
12 
    private $columns;
13 
 
14 
    private $rDiagonal = [];
15 
 
16 
    public function __construct(Matrix $matrix)
17 
    {
18 
        $this->qrMatrix = $matrix->toArray();
19 
        $this->rows = $matrix->rows;
20 
        $this->columns = $matrix->columns;
21 
 
22 
        $this->decompose();
23 
    }
24 
 
25 
    public function getHouseholdVectors(): Matrix
26 
    {
27 
        $householdVectors = [];
28 
        for ($row = 0; $row < $this->rows; ++$row) {
29 
            for ($column = 0; $column < $this->columns; ++$column) {
30 
                if ($row >= $column) {
31 
                    $householdVectors[$row][$column] = $this->qrMatrix[$row][$column];
32 
                } else {
33 
                    $householdVectors[$row][$column] = 0.0;
34 
                }
35 
            }
36 
        }
37 
 
38 
        return new Matrix($householdVectors);
39 
    }
40 
 
41 
    public function getQ(): Matrix
42 
    {
43 
        $qGrid = [];
44 
 
45 
        $rowCount = $this->rows;
46 
        for ($k = $this->columns - 1; $k >= 0; --$k) {
47 
            for ($i = 0; $i < $this->rows; ++$i) {
48 
                $qGrid[$i][$k] = 0.0;
49 
            }
50 
            $qGrid[$k][$k] = 1.0;
51 
            if ($this->columns > $this->rows) {
52 
                $qGrid = array_slice($qGrid, 0, $this->rows);
53 
            }
54 
 
55 
            for ($j = $k; $j < $this->columns; ++$j) {
56 
                if (isset($this->qrMatrix[$k], $this->qrMatrix[$k][$k]) && $this->qrMatrix[$k][$k] != 0.0) {
57 
                    $s = 0.0;
58 
                    for ($i = $k; $i < $this->rows; ++$i) {
59 
                        $s += $this->qrMatrix[$i][$k] * $qGrid[$i][$j];
60 
                    }
61 
                    $s = -$s / $this->qrMatrix[$k][$k];
62 
                    for ($i = $k; $i < $this->rows; ++$i) {
63 
                        $qGrid[$i][$j] += $s * $this->qrMatrix[$i][$k];
64 
                    }
65 
                }
66 
            }
67 
        }
68 
 
69 
        array_walk(
70 
            $qGrid,
71 
            function (&$row) use ($rowCount) {
72 
                $row = array_reverse($row);
73 
                $row = array_slice($row, 0, $rowCount);
74 
            }
75 
        );
76 
 
77 
        return new Matrix($qGrid);
78 
    }
79 
 
80 
    public function getR(): Matrix
81 
    {
82 
        $rGrid = [];
83 
 
84 
        for ($row = 0; $row < $this->columns; ++$row) {
85 
            for ($column = 0; $column < $this->columns; ++$column) {
86 
                if ($row < $column) {
87 
                    $rGrid[$row][$column] = $this->qrMatrix[$row][$column] ?? 0.0;
88 
                } elseif ($row === $column) {
89 
                    $rGrid[$row][$column] = $this->rDiagonal[$row] ?? 0.0;
90 
                } else {
91 
                    $rGrid[$row][$column] = 0.0;
92 
                }
93 
            }
94 
        }
95 
 
96 
        if ($this->columns > $this->rows) {
97 
            $rGrid = array_slice($rGrid, 0, $this->rows);
98 
        }
99 
 
100 
        return new Matrix($rGrid);
101 
    }
102 
 
103 
    private function hypo($a, $b): float
104 
    {
105 
        if (abs($a) > abs($b)) {
106 
            $r = $b / $a;
107 
            $r = abs($a) * sqrt(1 + $r * $r);
108 
        } elseif ($b != 0.0) {
109 
            $r = $a / $b;
110 
            $r = abs($b) * sqrt(1 + $r * $r);
111 
        } else {
112 
            $r = 0.0;
113 
        }
114 
 
115 
        return $r;
116 
    }
117 
 
118 
    /**
119 
     * QR Decomposition computed by Householder reflections.
120 
     */
121 
    private function decompose(): void
122 
    {
123 
        for ($k = 0; $k < $this->columns; ++$k) {
124 
            // Compute 2-norm of k-th column without under/overflow.
125 
            $norm = 0.0;
126 
            for ($i = $k; $i < $this->rows; ++$i) {
127 
                $norm = $this->hypo($norm, $this->qrMatrix[$i][$k]);
128 
            }
129 
            if ($norm != 0.0) {
130 
                // Form k-th Householder vector.
131 
                if ($this->qrMatrix[$k][$k] < 0.0) {
132 
                    $norm = -$norm;
133 
                }
134 
                for ($i = $k; $i < $this->rows; ++$i) {
135 
                    $this->qrMatrix[$i][$k] /= $norm;
136 
                }
137 
                $this->qrMatrix[$k][$k] += 1.0;
138 
                // Apply transformation to remaining columns.
139 
                for ($j = $k + 1; $j < $this->columns; ++$j) {
140 
                    $s = 0.0;
141 
                    for ($i = $k; $i < $this->rows; ++$i) {
142 
                        $s += $this->qrMatrix[$i][$k] * $this->qrMatrix[$i][$j];
143 
                    }
144 
                    $s = -$s / $this->qrMatrix[$k][$k];
145 
                    for ($i = $k; $i < $this->rows; ++$i) {
146 
                        $this->qrMatrix[$i][$j] += $s * $this->qrMatrix[$i][$k];
147 
                    }
148 
                }
149 
            }
150 
            $this->rDiagonal[$k] = -$norm;
151 
        }
152 
    }
153 
 
154 
    public function isFullRank(): bool
155 
    {
156 
        for ($j = 0; $j < $this->columns; ++$j) {
157 
            if ($this->rDiagonal[$j] == 0.0) {
158 
                return false;
159 
            }
160 
        }
161 
 
162 
        return true;
163 
    }
164 
 
165 
    /**
166 
     * Least squares solution of A*X = B.
167 
     *
168 
     * @param Matrix $B a Matrix with as many rows as A and any number of columns
169 
     *
170 
     * @throws Exception
171 
     *
172 
     * @return Matrix matrix that minimizes the two norm of Q*R*X-B
173 
     */
174 
    public function solve(Matrix $B): Matrix
175 
    {
176 
        if ($B->rows !== $this->rows) {
177 
            throw new Exception('Matrix row dimensions are not equal');
178 
        }
179 
 
180 
        if (!$this->isFullRank()) {
181 
            throw new Exception('Can only perform this operation on a full-rank matrix');
182 
        }
183 
 
184 
        // Compute Y = transpose(Q)*B
185 
        $Y = $this->getQ()->transpose()
186 
            ->multiply($B);
187 
        // Solve R*X = Y;
188 
        return $this->getR()->inverse()
189 
            ->multiply($Y);
190 
    }
191 
}