Chemical Data Processing Library C++ API - Version 1.1.1
SVDecomposition.hpp
Go to the documentation of this file.
1 /*
2  * SVDecomposition.hpp
3  *
4  * Copyright (C) 2003 Thomas Seidel <thomas.seidel@univie.ac.at>
5  *
6  * This library is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Lesser General Public
8  * License as published by the Free Software Foundation; either
9  * version 2 of the License, or (at your option) any later version.
10  *
11  * This library is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public License
17  * along with this library; see the file COPYING. If not, write to
18  * the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
19  * Boston, MA 02111-1307, USA.
20  */
21 
27 #ifndef CDPL_MATH_SVDECOMPOSITION_HPP
28 #define CDPL_MATH_SVDECOMPOSITION_HPP
29 
30 #include <cstddef>
31 #include <cmath>
32 #include <limits>
33 
34 #include "CDPL/Math/Check.hpp"
35 #include "CDPL/Math/CommonType.hpp"
36 #include "CDPL/Math/TypeTraits.hpp"
37 #include "CDPL/Math/MatrixProxy.hpp"
38 #include "CDPL/Math/Vector.hpp"
39 #include "CDPL/Math/SpecialFunctions.hpp"
40 #include "CDPL/Base/Exceptions.hpp"
41 
42 
43 namespace CDPL
44 {
45 
46  namespace Math
47  {
48 
69  template <typename A, typename W, typename V>
70  bool svDecompose(MatrixExpression<A>& a, VectorExpression<W>& w, MatrixExpression<V>& v, std::size_t max_iter = 0)
71  {
72  typedef typename CommonType<typename CommonType<typename A::ValueType, typename W::ValueType>::Type,
73  typename V::ValueType>::Type ValueType;
74 
75  typedef typename CommonType<typename CommonType<typename A::SizeType, typename W::SizeType>::Type,
76  typename V::SizeType>::Type SizeType;
77 
78  SizeType m = a().getSize1();
79  SizeType n = a().getSize2();
80 
81  CDPL_MATH_CHECK(SizeType(w().getSize()) >= n && SizeType(v().getSize1()) >= n && SizeType(v().getSize2()) >= n,
82  "svDecompose: Preconditions violated", Base::SizeError);
83 
84  Vector<ValueType> rv1(n);
85  SizeType l;
86  ValueType g, scale, anorm, s;
87 
88  g = scale = anorm = ValueType(0);
89 
90  for (SizeType i = 0; i < n; i++) {
91  l = i + 2;
92  rv1(i) = scale * g;
93  g = s = scale = ValueType(0);
94 
95  if (i < m) {
96  for (SizeType k = i; k < m; k++)
97  scale += TypeTraits<ValueType>::abs(a()(k, i));
98 
99  if (scale != ValueType(0)) {
100  for (SizeType k = i; k < m; k++) {
101  a()(k, i) /= scale;
102  s += a()(k, i) * a()(k, i);
103  }
104 
105  ValueType f = a()(i, i);
106 
107  g = -sign(TypeTraits<ValueType>::sqrt(s), f);
108 
109  ValueType h = f * g - s;
110 
111  a()(i, i) = f - g;
112 
113  for (SizeType j = l - 1; j < n; j++) {
114  s = ValueType(0);
115 
116  for (SizeType k = i; k < m; k++)
117  s += a()(k, i) * a()(k, j);
118 
119  f = s / h;
120 
121  for (SizeType k = i; k < m; k++)
122  a()(k, j) += f * a()(k, i);
123  }
124 
125  for (SizeType k = i; k < m; k++)
126  a()(k, i) *= scale;
127  }
128  }
129 
130  w()(i) = scale * g;
131  g = s = scale = ValueType(0);
132 
133  if (i + 1 <= m && i + 1 != n) {
134  for (SizeType k = l - 1; k < n; k++)
135  scale += TypeTraits<ValueType>::abs(a()(i, k));
136 
137  if (scale != ValueType(0)) {
138  for (SizeType k = l - 1; k < n; k++) {
139  a()(i, k) /= scale;
140  s += a()(i, k) * a()(i, k);
141  }
142 
143  ValueType f = a()(i, l - 1);
144 
145  g = -sign(TypeTraits<ValueType>::sqrt(s), f);
146 
147  ValueType h = f * g - s;
148 
149  a()(i, l - 1) = f - g;
150 
151  for (SizeType k = l - 1; k < n; k++)
152  rv1(k) = a()(i, k) / h;
153 
154  for (SizeType j = l - 1; j < m; j++) {
155  s = ValueType(0);
156 
157  for (SizeType k = l - 1; k < n; k++)
158  s += a()(j, k) * a()(i, k);
159 
160  for (SizeType k = l - 1; k < n; k++)
161  a()(j, k) += s * rv1(k);
162  }
163 
164  for (SizeType k = l - 1; k < n; k++)
165  a()(i, k) *= scale;
166  }
167  }
168 
169  anorm = std::max(anorm, (TypeTraits<ValueType>::abs(w()(i)) + TypeTraits<ValueType>::abs(rv1(i))));
170  }
171 
172  for (SizeType ic = n; ic > 0; ic--) {
173  SizeType i = ic - 1;
174 
175  if (i < n - 1) {
176  if (g != ValueType(0)) {
177  for (SizeType j = l; j < n; j++)
178  v()(j, i) = (a()(i, j) / a()(i, l)) / g;
179 
180  for (SizeType j = l; j < n; j++) {
181  s = ValueType(0);
182 
183  for (SizeType k = l; k < n; k++)
184  s += a()(i, k) * v()(k, j);
185 
186  for (SizeType k = l; k < n; k++)
187  v()(k, j) += s * v()(k, i);
188  }
189  }
190 
191  for (SizeType j = l; j < n; j++)
192  v()(i, j) = v()(j, i) = ValueType(0);
193  }
194 
195  v()(i, i) = ValueType(1);
196  g = rv1(i);
197  l = i;
198  }
199 
200  for (SizeType ic = std::min(m, n); ic > 0; ic--) {
201  SizeType i = ic - 1;
202 
203  l = ic;
204  g = w()(i);
205 
206  for (SizeType j = l; j < n; j++)
207  a()(i, j) = ValueType(0);
208 
209  if (g != ValueType(0)) {
210  g = ValueType(1) / g;
211 
212  for (SizeType j = l; j < n; j++) {
213  s = ValueType(0);
214 
215  for (SizeType k = l; k < m; k++)
216  s += a()(k, i) * a()(k, j);
217 
218  ValueType f = (s / a()(i, i)) * g;
219 
220  for (SizeType k = i; k < m; k++)
221  a()(k, j) += f * a()(k, i);
222  }
223 
224  for (SizeType j = i; j < m; j++)
225  a()(j, i) *= g;
226 
227  } else
228  for (SizeType j = i; j < m; j++)
229  a()(j, i) = ValueType(0);
230 
231  ++a()(i, i);
232  }
233 
234  for (SizeType kc = n; kc > 0; kc--) {
235  SizeType k = kc - 1;
236  SizeType nm = 0;
237 
238  for (std::size_t its = 0; max_iter == 0 || its < max_iter; its++) {
239  bool flag = true;
240 
241  for (SizeType lc = k + 1; lc > 0; lc--) {
242  l = lc - 1;
243  nm = l - 1;
244 
245  if (TypeTraits<ValueType>::abs(rv1(l)) + anorm == anorm) {
246  flag = false;
247  break;
248  }
249 
250  if (TypeTraits<ValueType>::abs(w()(nm)) + anorm == anorm)
251  break;
252  }
253 
254  if (flag) {
255  ValueType c = ValueType(0);
256 
257  s = ValueType(1);
258 
259  for (SizeType i = l; i < k + 1; i++) {
260  ValueType f = s * rv1(i);
261 
262  rv1(i) = c * rv1(i);
263 
264  if (TypeTraits<ValueType>::abs(f) + anorm == anorm)
265  break;
266 
267  g = w()(i);
268 
269  ValueType h = pythag(f, g);
270 
271  w()(i) = h;
272  h = ValueType(1) / h;
273  c = g * h;
274  s = -f * h;
275 
276  for (SizeType j = 0; j < m; j++) {
277  ValueType y = a()(j, nm);
278  ValueType z = a()(j, i);
279 
280  a()(j, nm) = y * c + z * s;
281  a()(j, i) = z * c - y * s;
282  }
283  }
284  }
285 
286  ValueType z = w()(k);
287 
288  if (l == k) {
289  if (z < ValueType(0)) {
290  w()(k) = -z;
291 
292  for (SizeType j = 0; j < n; j++)
293  v()(j, k) = -v()(j, k);
294  }
295  break;
296  }
297 
298  if (max_iter > 0 && its >= max_iter - 1)
299  return false;
300 
301  ValueType x = w()(l);
302  ValueType c = ValueType(1);
303 
304  nm = k - 1;
305 
306  ValueType y = w()(nm);
307 
308  g = rv1(nm);
309 
310  ValueType h = rv1(k);
311  ValueType f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
312 
313  g = pythag(f, ValueType(1));
314  f = ((x - z) * (x + z) + h * ((y / (f + sign(g, f))) - h)) / x;
315  s = ValueType(1);
316 
317  for (SizeType j = l; j <= nm; j++) {
318  SizeType i = j + 1;
319 
320  g = rv1(i);
321  y = w()(i);
322  h = s * g;
323  g = c * g;
324  z = pythag(f, h);
325  rv1(j) = z;
326  c = f / z;
327  s = h / z;
328  f = x * c + g * s;
329  g = g * c - x * s;
330  h = y * s;
331  y *= c;
332 
333  for (SizeType jj = 0; jj < n; jj++) {
334  x = v()(jj, j);
335  z = v()(jj, i);
336  v()(jj, j) = x * c + z * s;
337  v()(jj, i) = z * c - x * s;
338  }
339 
340  z = pythag(f, h);
341  w()(j) = z;
342 
343  if (z) {
344  z = ValueType(1) / z;
345  c = f * z;
346  s = h * z;
347  }
348 
349  f = c * g + s * y;
350  x = c * y - s * g;
351 
352  for (SizeType jj = 0; jj < m; jj++) {
353  y = a()(jj, j);
354  z = a()(jj, i);
355  a()(jj, j) = y * c + z * s;
356  a()(jj, i) = z * c - y * s;
357  }
358  }
359 
360  rv1(l) = ValueType(0);
361  rv1(k) = f;
362  w()(k) = x;
363  }
364  }
365 
366  // reordering
367 
368  SizeType inc = 1;
369  Vector<ValueType> su(m);
370  Vector<ValueType>& sv = rv1;
371 
372  do {
373  inc *= 3;
374  inc++;
375 
376  } while (inc <= n);
377 
378  do {
379  inc /= 3;
380 
381  for (SizeType i = inc; i < n; i++) {
382  ValueType sw = w()(i);
383 
384  for (SizeType k = 0; k < m; k++)
385  su(k) = a()(k, i);
386 
387  for (SizeType k = 0; k < n; k++)
388  sv(k) = v()(k, i);
389 
390  SizeType j = i;
391 
392  while (w()(j - inc) < sw) {
393  w()(j) = w()(j - inc);
394 
395  for (SizeType k = 0; k < m; k++)
396  a()(k, j) = a()(k, j - inc);
397 
398  for (SizeType k = 0; k < n; k++)
399  v()(k, j) = v()(k, j - inc);
400 
401  j -= inc;
402 
403  if (j < inc)
404  break;
405  }
406 
407  w()(j) = sw;
408 
409  for (SizeType k = 0; k < m; k++)
410  a()(k, j) = su(k);
411 
412  for (SizeType k = 0; k < n; k++)
413  v()(k, j) = sv(k);
414  }
415 
416  } while (inc > 1);
417 
418  for (SizeType k = 0; k < n; k++) {
419  SizeType s = 0;
420 
421  for (SizeType i = 0; i < m; i++)
422  if (a()(i, k) < typename A::ValueType(0))
423  s++;
424 
425  for (SizeType j = 0; j < n; j++)
426  if (v()(j, k) < typename V::ValueType(0))
427  s++;
428 
429  if (s > (m + n) / 2) {
430  for (SizeType i = 0; i < m; i++)
431  a()(i, k) = -a()(i, k);
432 
433  for (SizeType j = 0; j < n; j++)
434  v()(j, k) = -v()(j, k);
435  }
436  }
437 
438  return true;
439  }
440 
465  template <typename U, typename W, typename V, typename B, typename X>
466  void svSubstitute(const MatrixExpression<U>& u, const VectorExpression<W>& w, const MatrixExpression<V>& v,
467  const VectorExpression<B>& b, VectorExpression<X>& x)
468  {
469  typedef typename CommonType<typename CommonType<typename CommonType<typename CommonType<typename U::ValueType,
470  typename W::ValueType>::Type,
471  typename V::ValueType>::Type,
472  typename B::ValueType>::Type,
473  typename X::ValueType>::Type ValueType;
474 
475  typedef typename CommonType<typename U::SizeType, typename W::SizeType>::Type SizeType;
476 
477  SizeType m = u().getSize1();
478  SizeType n = u().getSize2();
479 
480  CDPL_MATH_CHECK(w().getSize() == n && v().getSize1() == n && v().getSize2() == n &&
481  b().getSize() == m && x().getSize() == n,
482  "svSubstitute: Preconditions violated", Base::SizeError);
483 
484  Vector<ValueType> tmp(n);
485  ValueType thresh = 0.5 * std::sqrt(m + n + 1.0) * w()(0) * std::numeric_limits<ValueType>::epsilon();
486 
487  for (SizeType i = 0; i < n; i++) { // Calculate trans(U) * B
488  ValueType s = ValueType(0);
489 
490  if (w()(i) > thresh) // Nonzero result only if w(i) > threshold
491  s = innerProd(column(u, i), b) / w()(i);
492 
493  tmp(i) = s;
494  }
495 
496  x().assign(prod(v, tmp));
497  }
498 
524  template <typename U, typename W, typename V, typename B, typename X>
525  void svSubstitute(const MatrixExpression<U>& u, const VectorExpression<W>& w, const MatrixExpression<V>& v,
526  const MatrixExpression<B>& b, MatrixExpression<X>& x)
527  {
528 
529  typedef typename CommonType<typename CommonType<typename CommonType<typename CommonType<typename U::ValueType,
530  typename W::ValueType>::Type,
531  typename V::ValueType>::Type,
532  typename B::ValueType>::Type,
533  typename X::ValueType>::Type ValueType;
534 
535  typedef typename CommonType<typename CommonType<typename CommonType<typename U::SizeType,
536  typename W::SizeType>::Type,
537  typename B::SizeType>::Type,
538  typename X::SizeType>::Type SizeType;
539 
540  SizeType m = u().getSize1();
541  SizeType n = u().getSize2();
542  SizeType p = b().getSize2();
543 
544  CDPL_MATH_CHECK(SizeType(w().getSize()) == n && SizeType(v().getSize1()) == n && SizeType(v().getSize2()) == n &&
545  SizeType(b().getSize1()) == m && SizeType(x().getSize1()) == n && SizeType(x().getSize2()) == p,
546  "svSubstitute: Preconditions violated", Base::SizeError);
547 
548  Vector<ValueType> tmp(n);
549  ValueType thresh = 0.5 * std::sqrt(m + n + 1.0) * w()(0) * std::numeric_limits<ValueType>::epsilon();
550 
551  for (SizeType i = 0; i < p; i++) { // For each column of B solve for the respective column of X
552  MatrixColumn<const B> b_col(column(b, i));
553 
554  for (SizeType j = 0; j < n; j++) { // Calculate trans(U) * B
555  ValueType s = ValueType(0);
556 
557  if (w()(j) > thresh) // Nonzero result only if w(j) > threshold
558  s = innerProd(column(u, j), b_col) / w()(j);
559 
560  tmp(j) = s;
561  }
562 
563  column(x, i).assign(prod(v, tmp));
564  }
565  }
566  } // namespace Math
567 } // namespace CDPL
568 
569 #endif // CDPL_MATH_SVDECOMPOSITION_HPP
CDPL::Chem::CIPDescriptor::p
const unsigned int p
Specifies that the stereocenter has p configuration.
Definition: CIPDescriptor.hpp:121
CDPL::Math::TypeTraits
Definition: TypeTraits.hpp:171
CDPL::Math::Vector< ValueType >
CDPL::Math::CommonType::Type
std::common_type< T1, T2 >::type Type
Definition: CommonType.hpp:43
CommonType.hpp
Common type deduction.
CDPL::Math::column
MatrixColumn< M > column(MatrixExpression< M > &e, typename MatrixColumn< M >::SizeType j)
Definition: MatrixProxy.hpp:730
CDPL::Math::VectorExpression
Definition: Expression.hpp:54
CDPL::Math::pythag
T pythag(const T &a, const T &b)
Computes without destructive underflow or overflow.
CDPL::Math::innerProd
VectorInnerProduct< E1, E2 >::ResultType innerProd(const VectorExpression< E1 > &e1, const VectorExpression< E2 > &e2)
Definition: VectorExpression.hpp:504
CDPL::Math::MatrixExpression
Definition: Expression.hpp:76
CDPL_MATH_CHECK
#define CDPL_MATH_CHECK(expr, msg, e)
Definition: Check.hpp:36
CDPL::Math::CommonType
Definition: CommonType.hpp:41
CDPL::Base::SizeError
Thrown to indicate that the size of a (multidimensional) array is not correct.
Definition: Base/Exceptions.hpp:133
CDPL::Math::MatrixColumn
Definition: MatrixProxy.hpp:196
TypeTraits.hpp
Definition of type traits.
CDPL::Math::prod
Matrix1VectorBinaryTraits< E1, E2, MatrixVectorProduct< E1, E2 > >::ResultType prod(const MatrixExpression< E1 > &e1, const VectorExpression< E2 > &e2)
Definition: MatrixExpression.hpp:833
CDPL::Chem::CIPDescriptor::m
const unsigned int m
Specifies that the stereocenter has m configuration.
Definition: CIPDescriptor.hpp:116
Exceptions.hpp
Definition of exception classes.
CDPL::Chem::CIPDescriptor::s
const unsigned int s
Specifies that the stereocenter has s configuration.
Definition: CIPDescriptor.hpp:81
SpecialFunctions.hpp
Provides miscellaneous special mathematical functions.
CDPL
The namespace of the Chemical Data Processing Library.
MatrixProxy.hpp
Definition of matrix proxy types.
CDPL::Math::sign
T1 sign(const T1 &a, const T2 &b)
Returns the magnitude of parameter a times the sign of parameter b.
Check.hpp
Definition of various preprocessor macros for error checking.
CDPL::Math::svSubstitute
void svSubstitute(const MatrixExpression< U > &u, const VectorExpression< W > &w, const MatrixExpression< V > &v, const VectorExpression< B > &b, VectorExpression< X > &x)
Solves for a vector where is given by its Singular Value Decomposition [WSVD].
Definition: SVDecomposition.hpp:466
CDPL::Math::svDecompose
bool svDecompose(MatrixExpression< A > &a, VectorExpression< W > &w, MatrixExpression< V > &v, std::size_t max_iter=0)
Computes the Singular Value Decomposition [WSVD] of a -dimensional matrix a.
Definition: SVDecomposition.hpp:70
Vector.hpp
Definition of vector data types.
Generated by doxygen 1.8.20