1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
|
---|
5 | //
|
---|
6 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
9 |
|
---|
10 | #ifndef EIGEN2_SVD_H
|
---|
11 | #define EIGEN2_SVD_H
|
---|
12 |
|
---|
13 | namespace Eigen {
|
---|
14 |
|
---|
15 | /** \ingroup SVD_Module
|
---|
16 | * \nonstableyet
|
---|
17 | *
|
---|
18 | * \class SVD
|
---|
19 | *
|
---|
20 | * \brief Standard SVD decomposition of a matrix and associated features
|
---|
21 | *
|
---|
22 | * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
|
---|
23 | *
|
---|
24 | * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N
|
---|
25 | * with \c M \>= \c N.
|
---|
26 | *
|
---|
27 | *
|
---|
28 | * \sa MatrixBase::SVD()
|
---|
29 | */
|
---|
30 | template<typename MatrixType> class SVD
|
---|
31 | {
|
---|
32 | private:
|
---|
33 | typedef typename MatrixType::Scalar Scalar;
|
---|
34 | typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
|
---|
35 |
|
---|
36 | enum {
|
---|
37 | PacketSize = internal::packet_traits<Scalar>::size,
|
---|
38 | AlignmentMask = int(PacketSize)-1,
|
---|
39 | MinSize = EIGEN_SIZE_MIN_PREFER_DYNAMIC(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime)
|
---|
40 | };
|
---|
41 |
|
---|
42 | typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVector;
|
---|
43 | typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> RowVector;
|
---|
44 |
|
---|
45 | typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MinSize> MatrixUType;
|
---|
46 | typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixVType;
|
---|
47 | typedef Matrix<Scalar, MinSize, 1> SingularValuesType;
|
---|
48 |
|
---|
49 | public:
|
---|
50 |
|
---|
51 | SVD() {} // a user who relied on compiler-generated default compiler reported problems with MSVC in 2.0.7
|
---|
52 |
|
---|
53 | SVD(const MatrixType& matrix)
|
---|
54 | : m_matU(matrix.rows(), (std::min)(matrix.rows(), matrix.cols())),
|
---|
55 | m_matV(matrix.cols(),matrix.cols()),
|
---|
56 | m_sigma((std::min)(matrix.rows(),matrix.cols()))
|
---|
57 | {
|
---|
58 | compute(matrix);
|
---|
59 | }
|
---|
60 |
|
---|
61 | template<typename OtherDerived, typename ResultType>
|
---|
62 | bool solve(const MatrixBase<OtherDerived> &b, ResultType* result) const;
|
---|
63 |
|
---|
64 | const MatrixUType& matrixU() const { return m_matU; }
|
---|
65 | const SingularValuesType& singularValues() const { return m_sigma; }
|
---|
66 | const MatrixVType& matrixV() const { return m_matV; }
|
---|
67 |
|
---|
68 | void compute(const MatrixType& matrix);
|
---|
69 | SVD& sort();
|
---|
70 |
|
---|
71 | template<typename UnitaryType, typename PositiveType>
|
---|
72 | void computeUnitaryPositive(UnitaryType *unitary, PositiveType *positive) const;
|
---|
73 | template<typename PositiveType, typename UnitaryType>
|
---|
74 | void computePositiveUnitary(PositiveType *positive, UnitaryType *unitary) const;
|
---|
75 | template<typename RotationType, typename ScalingType>
|
---|
76 | void computeRotationScaling(RotationType *unitary, ScalingType *positive) const;
|
---|
77 | template<typename ScalingType, typename RotationType>
|
---|
78 | void computeScalingRotation(ScalingType *positive, RotationType *unitary) const;
|
---|
79 |
|
---|
80 | protected:
|
---|
81 | /** \internal */
|
---|
82 | MatrixUType m_matU;
|
---|
83 | /** \internal */
|
---|
84 | MatrixVType m_matV;
|
---|
85 | /** \internal */
|
---|
86 | SingularValuesType m_sigma;
|
---|
87 | };
|
---|
88 |
|
---|
89 | /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix
|
---|
90 | *
|
---|
91 | * \note this code has been adapted from JAMA (public domain)
|
---|
92 | */
|
---|
93 | template<typename MatrixType>
|
---|
94 | void SVD<MatrixType>::compute(const MatrixType& matrix)
|
---|
95 | {
|
---|
96 | const int m = matrix.rows();
|
---|
97 | const int n = matrix.cols();
|
---|
98 | const int nu = (std::min)(m,n);
|
---|
99 | ei_assert(m>=n && "In Eigen 2.0, SVD only works for MxN matrices with M>=N. Sorry!");
|
---|
100 | ei_assert(m>1 && "In Eigen 2.0, SVD doesn't work on 1x1 matrices");
|
---|
101 |
|
---|
102 | m_matU.resize(m, nu);
|
---|
103 | m_matU.setZero();
|
---|
104 | m_sigma.resize((std::min)(m,n));
|
---|
105 | m_matV.resize(n,n);
|
---|
106 |
|
---|
107 | RowVector e(n);
|
---|
108 | ColVector work(m);
|
---|
109 | MatrixType matA(matrix);
|
---|
110 | const bool wantu = true;
|
---|
111 | const bool wantv = true;
|
---|
112 | int i=0, j=0, k=0;
|
---|
113 |
|
---|
114 | // Reduce A to bidiagonal form, storing the diagonal elements
|
---|
115 | // in s and the super-diagonal elements in e.
|
---|
116 | int nct = (std::min)(m-1,n);
|
---|
117 | int nrt = (std::max)(0,(std::min)(n-2,m));
|
---|
118 | for (k = 0; k < (std::max)(nct,nrt); ++k)
|
---|
119 | {
|
---|
120 | if (k < nct)
|
---|
121 | {
|
---|
122 | // Compute the transformation for the k-th column and
|
---|
123 | // place the k-th diagonal in m_sigma[k].
|
---|
124 | m_sigma[k] = matA.col(k).end(m-k).norm();
|
---|
125 | if (m_sigma[k] != 0.0) // FIXME
|
---|
126 | {
|
---|
127 | if (matA(k,k) < 0.0)
|
---|
128 | m_sigma[k] = -m_sigma[k];
|
---|
129 | matA.col(k).end(m-k) /= m_sigma[k];
|
---|
130 | matA(k,k) += 1.0;
|
---|
131 | }
|
---|
132 | m_sigma[k] = -m_sigma[k];
|
---|
133 | }
|
---|
134 |
|
---|
135 | for (j = k+1; j < n; ++j)
|
---|
136 | {
|
---|
137 | if ((k < nct) && (m_sigma[k] != 0.0))
|
---|
138 | {
|
---|
139 | // Apply the transformation.
|
---|
140 | Scalar t = matA.col(k).end(m-k).eigen2_dot(matA.col(j).end(m-k)); // FIXME dot product or cwise prod + .sum() ??
|
---|
141 | t = -t/matA(k,k);
|
---|
142 | matA.col(j).end(m-k) += t * matA.col(k).end(m-k);
|
---|
143 | }
|
---|
144 |
|
---|
145 | // Place the k-th row of A into e for the
|
---|
146 | // subsequent calculation of the row transformation.
|
---|
147 | e[j] = matA(k,j);
|
---|
148 | }
|
---|
149 |
|
---|
150 | // Place the transformation in U for subsequent back multiplication.
|
---|
151 | if (wantu & (k < nct))
|
---|
152 | m_matU.col(k).end(m-k) = matA.col(k).end(m-k);
|
---|
153 |
|
---|
154 | if (k < nrt)
|
---|
155 | {
|
---|
156 | // Compute the k-th row transformation and place the
|
---|
157 | // k-th super-diagonal in e[k].
|
---|
158 | e[k] = e.end(n-k-1).norm();
|
---|
159 | if (e[k] != 0.0)
|
---|
160 | {
|
---|
161 | if (e[k+1] < 0.0)
|
---|
162 | e[k] = -e[k];
|
---|
163 | e.end(n-k-1) /= e[k];
|
---|
164 | e[k+1] += 1.0;
|
---|
165 | }
|
---|
166 | e[k] = -e[k];
|
---|
167 | if ((k+1 < m) & (e[k] != 0.0))
|
---|
168 | {
|
---|
169 | // Apply the transformation.
|
---|
170 | work.end(m-k-1) = matA.corner(BottomRight,m-k-1,n-k-1) * e.end(n-k-1);
|
---|
171 | for (j = k+1; j < n; ++j)
|
---|
172 | matA.col(j).end(m-k-1) += (-e[j]/e[k+1]) * work.end(m-k-1);
|
---|
173 | }
|
---|
174 |
|
---|
175 | // Place the transformation in V for subsequent back multiplication.
|
---|
176 | if (wantv)
|
---|
177 | m_matV.col(k).end(n-k-1) = e.end(n-k-1);
|
---|
178 | }
|
---|
179 | }
|
---|
180 |
|
---|
181 |
|
---|
182 | // Set up the final bidiagonal matrix or order p.
|
---|
183 | int p = (std::min)(n,m+1);
|
---|
184 | if (nct < n)
|
---|
185 | m_sigma[nct] = matA(nct,nct);
|
---|
186 | if (m < p)
|
---|
187 | m_sigma[p-1] = 0.0;
|
---|
188 | if (nrt+1 < p)
|
---|
189 | e[nrt] = matA(nrt,p-1);
|
---|
190 | e[p-1] = 0.0;
|
---|
191 |
|
---|
192 | // If required, generate U.
|
---|
193 | if (wantu)
|
---|
194 | {
|
---|
195 | for (j = nct; j < nu; ++j)
|
---|
196 | {
|
---|
197 | m_matU.col(j).setZero();
|
---|
198 | m_matU(j,j) = 1.0;
|
---|
199 | }
|
---|
200 | for (k = nct-1; k >= 0; k--)
|
---|
201 | {
|
---|
202 | if (m_sigma[k] != 0.0)
|
---|
203 | {
|
---|
204 | for (j = k+1; j < nu; ++j)
|
---|
205 | {
|
---|
206 | Scalar t = m_matU.col(k).end(m-k).eigen2_dot(m_matU.col(j).end(m-k)); // FIXME is it really a dot product we want ?
|
---|
207 | t = -t/m_matU(k,k);
|
---|
208 | m_matU.col(j).end(m-k) += t * m_matU.col(k).end(m-k);
|
---|
209 | }
|
---|
210 | m_matU.col(k).end(m-k) = - m_matU.col(k).end(m-k);
|
---|
211 | m_matU(k,k) = Scalar(1) + m_matU(k,k);
|
---|
212 | if (k-1>0)
|
---|
213 | m_matU.col(k).start(k-1).setZero();
|
---|
214 | }
|
---|
215 | else
|
---|
216 | {
|
---|
217 | m_matU.col(k).setZero();
|
---|
218 | m_matU(k,k) = 1.0;
|
---|
219 | }
|
---|
220 | }
|
---|
221 | }
|
---|
222 |
|
---|
223 | // If required, generate V.
|
---|
224 | if (wantv)
|
---|
225 | {
|
---|
226 | for (k = n-1; k >= 0; k--)
|
---|
227 | {
|
---|
228 | if ((k < nrt) & (e[k] != 0.0))
|
---|
229 | {
|
---|
230 | for (j = k+1; j < nu; ++j)
|
---|
231 | {
|
---|
232 | Scalar t = m_matV.col(k).end(n-k-1).eigen2_dot(m_matV.col(j).end(n-k-1)); // FIXME is it really a dot product we want ?
|
---|
233 | t = -t/m_matV(k+1,k);
|
---|
234 | m_matV.col(j).end(n-k-1) += t * m_matV.col(k).end(n-k-1);
|
---|
235 | }
|
---|
236 | }
|
---|
237 | m_matV.col(k).setZero();
|
---|
238 | m_matV(k,k) = 1.0;
|
---|
239 | }
|
---|
240 | }
|
---|
241 |
|
---|
242 | // Main iteration loop for the singular values.
|
---|
243 | int pp = p-1;
|
---|
244 | int iter = 0;
|
---|
245 | Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
|
---|
246 | while (p > 0)
|
---|
247 | {
|
---|
248 | int k=0;
|
---|
249 | int kase=0;
|
---|
250 |
|
---|
251 | // Here is where a test for too many iterations would go.
|
---|
252 |
|
---|
253 | // This section of the program inspects for
|
---|
254 | // negligible elements in the s and e arrays. On
|
---|
255 | // completion the variables kase and k are set as follows.
|
---|
256 |
|
---|
257 | // kase = 1 if s(p) and e[k-1] are negligible and k<p
|
---|
258 | // kase = 2 if s(k) is negligible and k<p
|
---|
259 | // kase = 3 if e[k-1] is negligible, k<p, and
|
---|
260 | // s(k), ..., s(p) are not negligible (qr step).
|
---|
261 | // kase = 4 if e(p-1) is negligible (convergence).
|
---|
262 |
|
---|
263 | for (k = p-2; k >= -1; --k)
|
---|
264 | {
|
---|
265 | if (k == -1)
|
---|
266 | break;
|
---|
267 | if (ei_abs(e[k]) <= eps*(ei_abs(m_sigma[k]) + ei_abs(m_sigma[k+1])))
|
---|
268 | {
|
---|
269 | e[k] = 0.0;
|
---|
270 | break;
|
---|
271 | }
|
---|
272 | }
|
---|
273 | if (k == p-2)
|
---|
274 | {
|
---|
275 | kase = 4;
|
---|
276 | }
|
---|
277 | else
|
---|
278 | {
|
---|
279 | int ks;
|
---|
280 | for (ks = p-1; ks >= k; --ks)
|
---|
281 | {
|
---|
282 | if (ks == k)
|
---|
283 | break;
|
---|
284 | Scalar t = (ks != p ? ei_abs(e[ks]) : Scalar(0)) + (ks != k+1 ? ei_abs(e[ks-1]) : Scalar(0));
|
---|
285 | if (ei_abs(m_sigma[ks]) <= eps*t)
|
---|
286 | {
|
---|
287 | m_sigma[ks] = 0.0;
|
---|
288 | break;
|
---|
289 | }
|
---|
290 | }
|
---|
291 | if (ks == k)
|
---|
292 | {
|
---|
293 | kase = 3;
|
---|
294 | }
|
---|
295 | else if (ks == p-1)
|
---|
296 | {
|
---|
297 | kase = 1;
|
---|
298 | }
|
---|
299 | else
|
---|
300 | {
|
---|
301 | kase = 2;
|
---|
302 | k = ks;
|
---|
303 | }
|
---|
304 | }
|
---|
305 | ++k;
|
---|
306 |
|
---|
307 | // Perform the task indicated by kase.
|
---|
308 | switch (kase)
|
---|
309 | {
|
---|
310 |
|
---|
311 | // Deflate negligible s(p).
|
---|
312 | case 1:
|
---|
313 | {
|
---|
314 | Scalar f(e[p-2]);
|
---|
315 | e[p-2] = 0.0;
|
---|
316 | for (j = p-2; j >= k; --j)
|
---|
317 | {
|
---|
318 | Scalar t(numext::hypot(m_sigma[j],f));
|
---|
319 | Scalar cs(m_sigma[j]/t);
|
---|
320 | Scalar sn(f/t);
|
---|
321 | m_sigma[j] = t;
|
---|
322 | if (j != k)
|
---|
323 | {
|
---|
324 | f = -sn*e[j-1];
|
---|
325 | e[j-1] = cs*e[j-1];
|
---|
326 | }
|
---|
327 | if (wantv)
|
---|
328 | {
|
---|
329 | for (i = 0; i < n; ++i)
|
---|
330 | {
|
---|
331 | t = cs*m_matV(i,j) + sn*m_matV(i,p-1);
|
---|
332 | m_matV(i,p-1) = -sn*m_matV(i,j) + cs*m_matV(i,p-1);
|
---|
333 | m_matV(i,j) = t;
|
---|
334 | }
|
---|
335 | }
|
---|
336 | }
|
---|
337 | }
|
---|
338 | break;
|
---|
339 |
|
---|
340 | // Split at negligible s(k).
|
---|
341 | case 2:
|
---|
342 | {
|
---|
343 | Scalar f(e[k-1]);
|
---|
344 | e[k-1] = 0.0;
|
---|
345 | for (j = k; j < p; ++j)
|
---|
346 | {
|
---|
347 | Scalar t(numext::hypot(m_sigma[j],f));
|
---|
348 | Scalar cs( m_sigma[j]/t);
|
---|
349 | Scalar sn(f/t);
|
---|
350 | m_sigma[j] = t;
|
---|
351 | f = -sn*e[j];
|
---|
352 | e[j] = cs*e[j];
|
---|
353 | if (wantu)
|
---|
354 | {
|
---|
355 | for (i = 0; i < m; ++i)
|
---|
356 | {
|
---|
357 | t = cs*m_matU(i,j) + sn*m_matU(i,k-1);
|
---|
358 | m_matU(i,k-1) = -sn*m_matU(i,j) + cs*m_matU(i,k-1);
|
---|
359 | m_matU(i,j) = t;
|
---|
360 | }
|
---|
361 | }
|
---|
362 | }
|
---|
363 | }
|
---|
364 | break;
|
---|
365 |
|
---|
366 | // Perform one qr step.
|
---|
367 | case 3:
|
---|
368 | {
|
---|
369 | // Calculate the shift.
|
---|
370 | Scalar scale = (std::max)((std::max)((std::max)((std::max)(
|
---|
371 | ei_abs(m_sigma[p-1]),ei_abs(m_sigma[p-2])),ei_abs(e[p-2])),
|
---|
372 | ei_abs(m_sigma[k])),ei_abs(e[k]));
|
---|
373 | Scalar sp = m_sigma[p-1]/scale;
|
---|
374 | Scalar spm1 = m_sigma[p-2]/scale;
|
---|
375 | Scalar epm1 = e[p-2]/scale;
|
---|
376 | Scalar sk = m_sigma[k]/scale;
|
---|
377 | Scalar ek = e[k]/scale;
|
---|
378 | Scalar b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/Scalar(2);
|
---|
379 | Scalar c = (sp*epm1)*(sp*epm1);
|
---|
380 | Scalar shift(0);
|
---|
381 | if ((b != 0.0) || (c != 0.0))
|
---|
382 | {
|
---|
383 | shift = ei_sqrt(b*b + c);
|
---|
384 | if (b < 0.0)
|
---|
385 | shift = -shift;
|
---|
386 | shift = c/(b + shift);
|
---|
387 | }
|
---|
388 | Scalar f = (sk + sp)*(sk - sp) + shift;
|
---|
389 | Scalar g = sk*ek;
|
---|
390 |
|
---|
391 | // Chase zeros.
|
---|
392 |
|
---|
393 | for (j = k; j < p-1; ++j)
|
---|
394 | {
|
---|
395 | Scalar t = numext::hypot(f,g);
|
---|
396 | Scalar cs = f/t;
|
---|
397 | Scalar sn = g/t;
|
---|
398 | if (j != k)
|
---|
399 | e[j-1] = t;
|
---|
400 | f = cs*m_sigma[j] + sn*e[j];
|
---|
401 | e[j] = cs*e[j] - sn*m_sigma[j];
|
---|
402 | g = sn*m_sigma[j+1];
|
---|
403 | m_sigma[j+1] = cs*m_sigma[j+1];
|
---|
404 | if (wantv)
|
---|
405 | {
|
---|
406 | for (i = 0; i < n; ++i)
|
---|
407 | {
|
---|
408 | t = cs*m_matV(i,j) + sn*m_matV(i,j+1);
|
---|
409 | m_matV(i,j+1) = -sn*m_matV(i,j) + cs*m_matV(i,j+1);
|
---|
410 | m_matV(i,j) = t;
|
---|
411 | }
|
---|
412 | }
|
---|
413 | t = numext::hypot(f,g);
|
---|
414 | cs = f/t;
|
---|
415 | sn = g/t;
|
---|
416 | m_sigma[j] = t;
|
---|
417 | f = cs*e[j] + sn*m_sigma[j+1];
|
---|
418 | m_sigma[j+1] = -sn*e[j] + cs*m_sigma[j+1];
|
---|
419 | g = sn*e[j+1];
|
---|
420 | e[j+1] = cs*e[j+1];
|
---|
421 | if (wantu && (j < m-1))
|
---|
422 | {
|
---|
423 | for (i = 0; i < m; ++i)
|
---|
424 | {
|
---|
425 | t = cs*m_matU(i,j) + sn*m_matU(i,j+1);
|
---|
426 | m_matU(i,j+1) = -sn*m_matU(i,j) + cs*m_matU(i,j+1);
|
---|
427 | m_matU(i,j) = t;
|
---|
428 | }
|
---|
429 | }
|
---|
430 | }
|
---|
431 | e[p-2] = f;
|
---|
432 | iter = iter + 1;
|
---|
433 | }
|
---|
434 | break;
|
---|
435 |
|
---|
436 | // Convergence.
|
---|
437 | case 4:
|
---|
438 | {
|
---|
439 | // Make the singular values positive.
|
---|
440 | if (m_sigma[k] <= 0.0)
|
---|
441 | {
|
---|
442 | m_sigma[k] = m_sigma[k] < Scalar(0) ? -m_sigma[k] : Scalar(0);
|
---|
443 | if (wantv)
|
---|
444 | m_matV.col(k).start(pp+1) = -m_matV.col(k).start(pp+1);
|
---|
445 | }
|
---|
446 |
|
---|
447 | // Order the singular values.
|
---|
448 | while (k < pp)
|
---|
449 | {
|
---|
450 | if (m_sigma[k] >= m_sigma[k+1])
|
---|
451 | break;
|
---|
452 | Scalar t = m_sigma[k];
|
---|
453 | m_sigma[k] = m_sigma[k+1];
|
---|
454 | m_sigma[k+1] = t;
|
---|
455 | if (wantv && (k < n-1))
|
---|
456 | m_matV.col(k).swap(m_matV.col(k+1));
|
---|
457 | if (wantu && (k < m-1))
|
---|
458 | m_matU.col(k).swap(m_matU.col(k+1));
|
---|
459 | ++k;
|
---|
460 | }
|
---|
461 | iter = 0;
|
---|
462 | p--;
|
---|
463 | }
|
---|
464 | break;
|
---|
465 | } // end big switch
|
---|
466 | } // end iterations
|
---|
467 | }
|
---|
468 |
|
---|
469 | template<typename MatrixType>
|
---|
470 | SVD<MatrixType>& SVD<MatrixType>::sort()
|
---|
471 | {
|
---|
472 | int mu = m_matU.rows();
|
---|
473 | int mv = m_matV.rows();
|
---|
474 | int n = m_matU.cols();
|
---|
475 |
|
---|
476 | for (int i=0; i<n; ++i)
|
---|
477 | {
|
---|
478 | int k = i;
|
---|
479 | Scalar p = m_sigma.coeff(i);
|
---|
480 |
|
---|
481 | for (int j=i+1; j<n; ++j)
|
---|
482 | {
|
---|
483 | if (m_sigma.coeff(j) > p)
|
---|
484 | {
|
---|
485 | k = j;
|
---|
486 | p = m_sigma.coeff(j);
|
---|
487 | }
|
---|
488 | }
|
---|
489 | if (k != i)
|
---|
490 | {
|
---|
491 | m_sigma.coeffRef(k) = m_sigma.coeff(i); // i.e.
|
---|
492 | m_sigma.coeffRef(i) = p; // swaps the i-th and the k-th elements
|
---|
493 |
|
---|
494 | int j = mu;
|
---|
495 | for(int s=0; j!=0; ++s, --j)
|
---|
496 | std::swap(m_matU.coeffRef(s,i), m_matU.coeffRef(s,k));
|
---|
497 |
|
---|
498 | j = mv;
|
---|
499 | for (int s=0; j!=0; ++s, --j)
|
---|
500 | std::swap(m_matV.coeffRef(s,i), m_matV.coeffRef(s,k));
|
---|
501 | }
|
---|
502 | }
|
---|
503 | return *this;
|
---|
504 | }
|
---|
505 |
|
---|
506 | /** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A.
|
---|
507 | * The parts of the solution corresponding to zero singular values are ignored.
|
---|
508 | *
|
---|
509 | * \sa MatrixBase::svd(), LU::solve(), LLT::solve()
|
---|
510 | */
|
---|
511 | template<typename MatrixType>
|
---|
512 | template<typename OtherDerived, typename ResultType>
|
---|
513 | bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* result) const
|
---|
514 | {
|
---|
515 | ei_assert(b.rows() == m_matU.rows());
|
---|
516 |
|
---|
517 | Scalar maxVal = m_sigma.cwise().abs().maxCoeff();
|
---|
518 | for (int j=0; j<b.cols(); ++j)
|
---|
519 | {
|
---|
520 | Matrix<Scalar,MatrixUType::RowsAtCompileTime,1> aux = m_matU.transpose() * b.col(j);
|
---|
521 |
|
---|
522 | for (int i = 0; i <m_matU.cols(); ++i)
|
---|
523 | {
|
---|
524 | Scalar si = m_sigma.coeff(i);
|
---|
525 | if (ei_isMuchSmallerThan(ei_abs(si),maxVal))
|
---|
526 | aux.coeffRef(i) = 0;
|
---|
527 | else
|
---|
528 | aux.coeffRef(i) /= si;
|
---|
529 | }
|
---|
530 |
|
---|
531 | result->col(j) = m_matV * aux;
|
---|
532 | }
|
---|
533 | return true;
|
---|
534 | }
|
---|
535 |
|
---|
536 | /** Computes the polar decomposition of the matrix, as a product unitary x positive.
|
---|
537 | *
|
---|
538 | * If either pointer is zero, the corresponding computation is skipped.
|
---|
539 | *
|
---|
540 | * Only for square matrices.
|
---|
541 | *
|
---|
542 | * \sa computePositiveUnitary(), computeRotationScaling()
|
---|
543 | */
|
---|
544 | template<typename MatrixType>
|
---|
545 | template<typename UnitaryType, typename PositiveType>
|
---|
546 | void SVD<MatrixType>::computeUnitaryPositive(UnitaryType *unitary,
|
---|
547 | PositiveType *positive) const
|
---|
548 | {
|
---|
549 | ei_assert(m_matU.cols() == m_matV.cols() && "Polar decomposition is only for square matrices");
|
---|
550 | if(unitary) *unitary = m_matU * m_matV.adjoint();
|
---|
551 | if(positive) *positive = m_matV * m_sigma.asDiagonal() * m_matV.adjoint();
|
---|
552 | }
|
---|
553 |
|
---|
554 | /** Computes the polar decomposition of the matrix, as a product positive x unitary.
|
---|
555 | *
|
---|
556 | * If either pointer is zero, the corresponding computation is skipped.
|
---|
557 | *
|
---|
558 | * Only for square matrices.
|
---|
559 | *
|
---|
560 | * \sa computeUnitaryPositive(), computeRotationScaling()
|
---|
561 | */
|
---|
562 | template<typename MatrixType>
|
---|
563 | template<typename UnitaryType, typename PositiveType>
|
---|
564 | void SVD<MatrixType>::computePositiveUnitary(UnitaryType *positive,
|
---|
565 | PositiveType *unitary) const
|
---|
566 | {
|
---|
567 | ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
|
---|
568 | if(unitary) *unitary = m_matU * m_matV.adjoint();
|
---|
569 | if(positive) *positive = m_matU * m_sigma.asDiagonal() * m_matU.adjoint();
|
---|
570 | }
|
---|
571 |
|
---|
572 | /** decomposes the matrix as a product rotation x scaling, the scaling being
|
---|
573 | * not necessarily positive.
|
---|
574 | *
|
---|
575 | * If either pointer is zero, the corresponding computation is skipped.
|
---|
576 | *
|
---|
577 | * This method requires the Geometry module.
|
---|
578 | *
|
---|
579 | * \sa computeScalingRotation(), computeUnitaryPositive()
|
---|
580 | */
|
---|
581 | template<typename MatrixType>
|
---|
582 | template<typename RotationType, typename ScalingType>
|
---|
583 | void SVD<MatrixType>::computeRotationScaling(RotationType *rotation, ScalingType *scaling) const
|
---|
584 | {
|
---|
585 | ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
|
---|
586 | Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1
|
---|
587 | Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma);
|
---|
588 | sv.coeffRef(0) *= x;
|
---|
589 | if(scaling) scaling->lazyAssign(m_matV * sv.asDiagonal() * m_matV.adjoint());
|
---|
590 | if(rotation)
|
---|
591 | {
|
---|
592 | MatrixType m(m_matU);
|
---|
593 | m.col(0) /= x;
|
---|
594 | rotation->lazyAssign(m * m_matV.adjoint());
|
---|
595 | }
|
---|
596 | }
|
---|
597 |
|
---|
598 | /** decomposes the matrix as a product scaling x rotation, the scaling being
|
---|
599 | * not necessarily positive.
|
---|
600 | *
|
---|
601 | * If either pointer is zero, the corresponding computation is skipped.
|
---|
602 | *
|
---|
603 | * This method requires the Geometry module.
|
---|
604 | *
|
---|
605 | * \sa computeRotationScaling(), computeUnitaryPositive()
|
---|
606 | */
|
---|
607 | template<typename MatrixType>
|
---|
608 | template<typename ScalingType, typename RotationType>
|
---|
609 | void SVD<MatrixType>::computeScalingRotation(ScalingType *scaling, RotationType *rotation) const
|
---|
610 | {
|
---|
611 | ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
|
---|
612 | Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1
|
---|
613 | Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma);
|
---|
614 | sv.coeffRef(0) *= x;
|
---|
615 | if(scaling) scaling->lazyAssign(m_matU * sv.asDiagonal() * m_matU.adjoint());
|
---|
616 | if(rotation)
|
---|
617 | {
|
---|
618 | MatrixType m(m_matU);
|
---|
619 | m.col(0) /= x;
|
---|
620 | rotation->lazyAssign(m * m_matV.adjoint());
|
---|
621 | }
|
---|
622 | }
|
---|
623 |
|
---|
624 |
|
---|
625 | /** \svd_module
|
---|
626 | * \returns the SVD decomposition of \c *this
|
---|
627 | */
|
---|
628 | template<typename Derived>
|
---|
629 | inline SVD<typename MatrixBase<Derived>::PlainObject>
|
---|
630 | MatrixBase<Derived>::svd() const
|
---|
631 | {
|
---|
632 | return SVD<PlainObject>(derived());
|
---|
633 | }
|
---|
634 |
|
---|
635 | } // end namespace Eigen
|
---|
636 |
|
---|
637 | #endif // EIGEN2_SVD_H
|
---|