Bayes++ Bayesian Filtering Classes Release 2014.5 - Copyright (c) 2003,2004,2005,2006,2011,2012,2014 Michael Stevens
uBLASmatrix.hpp
Go to the documentation of this file.
1/*
2 * Bayes++ the Bayesian Filtering Library
3 * Copyright (c) 2002 Michael Stevens
4 * See accompanying Bayes++.htm for terms and conditions of use.
5 *
6 * $Id$
7 */
8
9/*
10 * Common type independent uBlas interface
11 * Should be include after base types have been defined
12 *
13 * Everything in namespace Bayes_filter_matrix is intended to support the matrix storage
14 * and algebra requirements of the library. Therefore the interfaces and implementation is
15 * not intended to be stable. Nor is this a general purpose adaptor for uBLAS
16 *
17 * Note on row_major matrices
18 * The implementation uses row major extensively. The computation of symmetric products P*P' is
19 * most efficient with row operations. These products are used extensively so the default
20 * is to use row_major matrices
21 */
22
23/* Filter Matrix Namespace */
25{
26 // Allow use a few functions in own namespace (particular useful for compilers with Konig lookup)
27using ublas::row;
28using ublas::column;
29using ublas::trans;
30using ublas::prod; // These do not apply to the templated prod<temp> funtions
31using ublas::inner_prod;
32using ublas::outer_prod;
33
34enum EmptyTag {Empty}; // Tag type used for empty matrix constructor
35
36 // Old compiler workaround removed from new uBLAS
37#ifndef BOOST_UBLAS_TYPENAME
38#define BOOST_UBLAS_TYPENAME typename
39#endif
40
41#if (BOOST_VERSION >= 103100)
42using ublas::noalias;
43#else
44namespace detail
45{
46 // Assignment proxy.
47 // Provides temporary free assignment when LHS has no alias on RHS
48 template<class C>
50 public:
51 BOOST_UBLAS_INLINE
52 noalias_proxy (C& lval):
53 lval_ (lval) {}
54
55 template <class E>
56 BOOST_UBLAS_INLINE
57 void operator= (const E& e) {
58 lval_.assign (e);
59 }
60
61 template <class E>
62 BOOST_UBLAS_INLINE
63 void operator+= (const E& e) {
64 lval_.plus_assign (e);
65 }
66
67 template <class E>
68 BOOST_UBLAS_INLINE
69 void operator-= (const E& e) {
70 lval_.minus_assign (e);
71 }
72
73 private: // nonassignable
74 void operator=( const noalias_proxy& );
75 private:
76 C& lval_;
77 };
78
79}//namespace detail
80
81// Improve syntax of efficient assignment where no aliases of LHS appear on the RHS
82// noalias(lhs) = rhs_expression
83template <class E>
84BOOST_UBLAS_INLINE
85detail::noalias_proxy<E> noalias (ublas::matrix_expression<E>& lvalue) {
86 return detail::noalias_proxy<E> (lvalue() );
87}
88template <class E>
89BOOST_UBLAS_INLINE
90detail::noalias_proxy<E> noalias (ublas::vector_expression<E>& lvalue) {
91 return detail::noalias_proxy<E> (lvalue() );
92}
93
94#endif
95
96
97namespace detail // Lots of implementation detail
98{
99
100/*
101 * Filter Vec type
102 */
103template <class VecBase>
104class FMVec : public VecBase
105{
106public:
107 typedef typename VecBase::value_type value_type;
108 typedef typename VecBase::vector_temporary_type vector_temporary_type;
109
110 // No Default Constructor. Empty creation is very error prone
111 explicit FMVec(EmptyTag) : VecBase()
112 {} // Empty constructor
113 explicit FMVec(std::size_t size) : VecBase(size)
114 {} // Normal sized constructor
115 FMVec(const FMVec& c) : VecBase(static_cast<const VecBase&>(c))
116 {} // Copy constructor
117 template <class E>
118 explicit FMVec(const ublas::vector_expression<E>& e) : VecBase(e)
119 {} // vector_expression copy constructor
120 template <class E>
121 FMVec(const ublas::matrix_column<E>& e) : VecBase(e)
122 {} // conversion copy constructor, hides the implicit copy required for matrix column access
123
124 template <class E>
125 FMVec& operator= (const ublas::vector_expression<E>& r)
126 { // Expression assignment; may be dependent on r
127 VecBase::operator=(r);
128 return *this;
129 }
130 FMVec& operator= (const FMVec& r)
131 { // Vector assignment; independent
132 VecBase::assign(r);
133 return *this;
134 }
135
136 // Sub-range selection operators
137 const ublas::vector_range<const VecBase> sub_range(std::size_t b, std::size_t e) const
138 {
139 return ublas::vector_range<const VecBase>(*this, ublas::range(b,e));
140 }
141 ublas::vector_range<VecBase> sub_range(std::size_t b, std::size_t e)
142 {
143 return ublas::vector_range<VecBase>(*this, ublas::range(b,e));
144 }
145};
146
147
148/*
149 * Filter Matrix class template. Augmentation for uBlas MatrixBase
150 */
151template <class MatrixBase>
152class FMMatrix : public MatrixBase
153{
154public:
155 typedef typename MatrixBase::value_type value_type;
156 typedef typename MatrixBase::vector_temporary_type vector_temporary_type;
157 typedef typename MatrixBase::matrix_temporary_type matrix_temporary_type;
158
159 // No Default Constructor. Empty creation is very error prone
160 explicit FMMatrix(EmptyTag) : MatrixBase()
161 {} // Empty constructor
162 FMMatrix(std::size_t size1, std::size_t size2) : MatrixBase(size1,size2)
163 {} // Normal sized constructor
164 FMMatrix(const FMMatrix& c) : MatrixBase(static_cast<const MatrixBase&>(c))
165 {} // Copy constructor
166 template <class E>
167 explicit FMMatrix(const ublas::matrix_expression<E>& e) : MatrixBase(e)
168 {} // matrix_expression copy constructor
169
170 template <class E>
171 FMMatrix& operator= (const ublas::matrix_expression<E>& r)
172 { // Expression assignment; may be dependent on r
173 MatrixBase::operator=(r);
174 return *this;
175 }
176 FMMatrix& operator= (const FMMatrix& r)
177 { // Matrix assignment; independent
178 MatrixBase::assign (r);
179 return *this;
180 }
181
182 // Row,Column vector proxies
183 typedef ublas::matrix_row<FMMatrix> Row;
184 typedef const ublas::matrix_row<const FMMatrix> const_Row;
185 typedef ublas::matrix_column<FMMatrix> Column;
186 typedef const ublas::matrix_column<const FMMatrix> const_Column;
187
188 // Vector proxies from iterators - static members dependent on MatrixBase type
189 // ri() returns container associated with iterator. static_cast required as typeof(ri()) may not be typeof(MM)
190 static Row rowi(const typename MatrixBase::iterator1& ri)
191 {
192 return Row(static_cast<FMMatrix&>(ri()), ri.index1());
193 }
194 static const_Row rowi(const typename MatrixBase::const_iterator1& ri)
195 {
196 return const_Row(static_cast<const FMMatrix&>(ri()), ri.index1());
197 }
198 static Column columni(const typename MatrixBase::iterator2& ci)
199 {
200 return Column(static_cast<FMMatrix&>(ci()), ci.index2());
201 }
202 static const_Column columni(const typename MatrixBase::const_iterator2& ci)
203 {
204 return const_Column(static_cast<const FMMatrix&>(ci()), ci.index2());
205 }
206
207 // Sub-range selection operators
208 ublas::matrix_range<const MatrixBase>
209 sub_matrix(std::size_t s1, std::size_t e1, std::size_t s2, std::size_t e2) const
210 {
211 return ublas::matrix_range<const MatrixBase> (*this, ublas::range(s1,e1), ublas::range(s2,e2));
212 }
213 ublas::matrix_range<MatrixBase>
214 sub_matrix(std::size_t s1, std::size_t e1, std::size_t s2, std::size_t e2)
215 {
216 return ublas::matrix_range<MatrixBase> (*this, ublas::range(s1,e1), ublas::range(s2,e2));
217 }
218
219 // Requires boost_1.30.0 which has a generalised matrix_vector_slice
220 ublas::matrix_vector_slice<const MatrixBase>
221 sub_column(std::size_t s1, std::size_t e1, std::size_t s2) const
222 // Column vector s2 with rows [s1,e1)
223 {
224 return ublas::matrix_vector_slice<const MatrixBase> (*this, ublas::slice(s1,1,e1-s1), ublas::slice(s2,0,e1-s1));
225 }
226 ublas::matrix_vector_slice<MatrixBase>
227 sub_column(std::size_t s1, std::size_t e1, std::size_t s2)
228 // Column vector s2 with rows [s1,e1)
229 {
230 return ublas::matrix_vector_slice<MatrixBase> (*this, ublas::slice(s1,1,e1-s1), ublas::slice(s2,0,e1-s1));
231 }
232};
233
234
235/*
236 * Helper template to allow member construction before base class
237 * Boost version does not work as it passes by value
238 */
239template <typename MemberType>
240class BaseFromMember
241{
242protected:
243 MemberType member;
244 explicit BaseFromMember() : member()
245 {}
246
247 template <typename T1>
248 explicit BaseFromMember( const T1& x1 ) : member( x1 )
249 {}
250
251 template <typename T1, typename T2>
252 explicit BaseFromMember( const T1& x1, const T2& x2 ) : member( x1, x2 )
253 {}
254};
255
256
257/*
258 * We require static type conversion between Symmetric matrices and equivalent row major matrices
259 * Therefore we create symmetric matrix types, using a MatrixBase for storage
260 * and wraps this in a symmetric_adaptor
261 */
262template <class MatrixBase>
263class SymMatrixWrapper :
264 private BaseFromMember<MatrixBase>, // allow construction of MatrixBase member before symmetric_adaptor
265 public ublas::symmetric_adaptor<MatrixBase, ublas::upper>
266{
267 typedef BaseFromMember<MatrixBase> matrix_type;
268 typedef ublas::symmetric_adaptor<MatrixBase, ublas::upper> symadaptor_type;
269public:
270 typedef typename MatrixBase::value_type value_type;
271 typedef typename MatrixBase::vector_temporary_type vector_temporary_type;
272 typedef typename MatrixBase::matrix_temporary_type matrix_temporary_type;
273
274 SymMatrixWrapper () : matrix_type(), symadaptor_type(matrix_type::member)
275 {}
276 SymMatrixWrapper (std::size_t size1, std::size_t size2) : matrix_type(size1,size2), symadaptor_type(matrix_type::member)
277 {} // Normal sized constructor
278 explicit SymMatrixWrapper (const SymMatrixWrapper& r) : matrix_type(reinterpret_cast<const MatrixBase&>(r)), symadaptor_type(matrix_type::member)
279 {} // Explicit copy construction referencing the copy reinterpreted as a MatrixBase
280 template <class E>
281 explicit SymMatrixWrapper (const ublas::matrix_expression<E>& e) : matrix_type(e), symadaptor_type(matrix_type::member)
282 {} // Explicit matrix_expression conversion constructor
283
284 template <class E>
285 SymMatrixWrapper& operator=(const ublas::matrix_expression<E>& r)
286 {
287 symadaptor_type::operator=(r);
288 return *this;
289 }
290
291 // Conversions straight to a FMMatrix, equivalent to a RowMatrix types
292 const FMMatrix<MatrixBase>& asRowMatrix() const
293 {
294 return static_cast<const FMMatrix<MatrixBase>& >(matrix_type::member);
295 }
296 FMMatrix<MatrixBase>& asRowMatrix()
297 {
298 return static_cast<FMMatrix<MatrixBase>& >(matrix_type::member);
299 }
300
301 // Matrix storage members
302 void clear()
303 { matrix_type::member.clear();
304 }
305 void resize(std::size_t nsize1, std::size_t nsize2, bool preserve = true)
306 {
307 matrix_type::member.resize(nsize1, nsize2, preserve);
308 }
309};
310
311}//namespace detail
312
313
314
315/*
316 * Vector / Matrix types
317 * Finally the definitions !
318 */
319using detail::FMVec; // Template class for template parameter matching
320using detail::FMMatrix;
321
322 // Default types
323typedef FMVec<detail::BaseVector> Vec;
324typedef FMMatrix<detail::BaseRowMatrix> RowMatrix;
326typedef FMMatrix<detail::BaseColMatrix> ColMatrix;
327typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseRowMatrix> > SymMatrix;
328typedef FMMatrix<detail::BaseUpperTriMatrix> UTriMatrix;
329typedef FMMatrix<detail::BaseLowerTriMatrix> LTriMatrix;
330typedef FMMatrix<detail::BaseDiagMatrix> DiagMatrix;
331
332 // Explicitly dense types
333typedef FMVec<detail::BaseDenseVector> DenseVec;
334typedef FMMatrix<detail::BaseDenseRowMatrix> DenseRowMatrix;
336typedef FMMatrix<detail::BaseDenseColMatrix> DenseColMatrix;
337typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseDenseRowMatrix> > DenseSymMatrix;
338typedef FMMatrix<detail::BaseDenseUpperTriMatrix> DenseUTriMatrix;
339typedef FMMatrix<detail::BaseDenseLowerTriMatrix> DenseLTriMatrix;
340typedef FMMatrix<detail::BaseDenseDiagMatrix> DenseDiagMatrix;
341
342 // Explicitly sparse types (any of the gappy types)
343#ifdef BAYES_FILTER_GAPPY
344typedef FMVec<detail::BaseSparseVector> SparseVec;
345typedef FMMatrix<detail::BaseDenseRowMatrix> SparseRowMatrix;
346typedef SparseRowMatrix SparseMatrix;
347typedef FMMatrix<detail::BaseSparseColMatrix> SparseColMatrix;
348typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseSparseRowMatrix> > SparseSymMatrix;
349#endif
350
351
352/*
353 * Matrix Adaptors, simply hide the uBLAS details
354 */
355template <class M>
356const ublas::triangular_adaptor<const M, ublas::upper>
357 UpperTri(const M& m)
358/*
359 * View Upper triangle of m
360 * ISSUE VC7 cannot cope with UTriMatrix::functor1_type
361 */
362{
363 return ublas::triangular_adaptor<const M, ublas::upper>(m);
364}
365
366template <class M>
367const ublas::triangular_adaptor<const M, ublas::lower>
368 LowerTri(const M& m)
369/*
370 * View Lower triangle of m
371 */
372{
373 return ublas::triangular_adaptor<const M, ublas::lower>(m);
374}
375
376
377
378/*
379 * Matrix Support Operations
380 */
381template <class Base>
382ublas::matrix_vector_range<FMMatrix<Base> >
383 diag(FMMatrix<Base>& M, std::size_t n)
384{ // Return a vector proxy to the first n diagonal elements of M
385 return ublas::matrix_vector_range<FMMatrix<Base> >(M, ublas::range(0,n), ublas::range(0,n));
386}
387
388template <class Base>
389const ublas::matrix_vector_range<const FMMatrix<Base> >
390 diag(const FMMatrix<Base>& M, std::size_t n)
391{ // Return a vector proxy to the first n diagonal elements of M
392 return ublas::matrix_vector_range<const FMMatrix<Base> >(M, ublas::range(0,n), ublas::range(0,n));
393}
394
395template <class Base>
396ublas::matrix_vector_range<FMMatrix<Base> >
397 diag(FMMatrix<Base>& M)
398{ // Return a vector proxy to the diagonal elements of M
399 const std::size_t common_size = std::min(M.size1(),M.size2());
400 return ublas::matrix_vector_range<FMMatrix<Base> >(M, ublas::range(0,common_size), ublas::range(0,common_size));
401}
402
403template <class Base>
404const ublas::matrix_vector_range<const FMMatrix<Base> >
405 diag(const FMMatrix<Base>& M)
406{ // Return a vector proxy to the diagonal elements of M
407 const std::size_t common_size = std::min(M.size1(),M.size2());
408 return ublas::matrix_vector_range<const FMMatrix<Base> >(M, ublas::range(0,common_size), ublas::range(0,common_size));
409}
410
411template <class Base>
412void identity(FMMatrix<Base>& I)
413{ // Set I to generalised Identity matrix. Clear I and set diag(I) to one
414 I.clear();
415 // Set common diagonal elements
416 std::size_t common_size = std::min(I.size1(),I.size2());
417 typedef typename Base::value_type Base_value_type;
418 diag(I).assign (ublas::scalar_vector<Base_value_type>(common_size, Base_value_type(1)));
419}
420
421
422
423/*
424 * Symmetric Positive (Semi) Definite multiplication: X*S*X' and X'*S*X
425 * The name is slightly misleading. The result is actually only PD if S is PD
426 * Algorithms are intended to exploit the symmetry of the result
427 * and also where possible the row by row multiplication inherent in X*X'
428 */
429
430namespace detail { // mult_SPD now an implementation detail
431
432template <class MatrixX>
433void mult_SPD (const MatrixX& X, const Vec& s, SymMatrix& P)
434/*
435 * Symmetric Positive (Semi) Definite multiply: P += X*diag_matrix(s)*X'
436 */
437{
438 Vec::const_iterator si, send = s.end();
439 typename MatrixX::const_iterator1 Xa = X.begin1();
440 const typename MatrixX::const_iterator1 Xend = X.end1();
441 typename MatrixX::const_iterator1 Xb;
442
443 // P(a,b) = X.row(a) * X.row(b)
444 for (; Xa != Xend; ++Xa) // Iterate Rows
445 {
446 typename MatrixX::const_Row Xav = MatrixX::rowi(Xa);
447 Xb = Xa; // Start at the row Xa only one triangle of symmetric result required
448 for (; Xb != Xend; ++Xb)
449 {
450 SymMatrix::value_type p = 0; // Triple vector inner product
451 typename MatrixX::const_Row Xbv = MatrixX::rowi(Xb);
452 for (si = s.begin(); si != send; ++si) {
453 Vec::size_type i = si.index();
454 p += Xav[i] * (*si) * Xbv[i];
455 }
456 P(Xa.index1(),Xb.index1()) += p;
457 }
458 }
459}
460
461template <class MatrixX>
462void mult_SPDT (const MatrixX& X, const Vec& s, SymMatrix& P)
463/*
464 * Symmetric Positive (Semi) Definite multiply: P += X'*diag_matrix(s)*X
465 */
466{
467 Vec::const_iterator si, send = s.end();
468 typename MatrixX::const_iterator2 Xa = X.begin2();
469 const typename MatrixX::const_iterator2 Xend = X.end2();
470 typename MatrixX::const_iterator2 Xb;
471
472 // P(a,b) = X.col(a) * X.col(b)
473 for (; Xa != Xend; ++Xa) // Iterate vectors
474 {
475 typename MatrixX::const_Column Xav = MatrixX::columni(Xa);
476 Xb = Xa; // Start at the row Xa only one triangle of symmetric result required
477 for (; Xb != Xend; ++Xb)
478 { // Triple vector inner product
479 SymMatrix::value_type p = 0;
480 typename MatrixX::const_Column Xbv = MatrixX::columni(Xb);
481 for (si = s.begin(); si != send; ++si) {
482 Vec::size_type i = si.index();
483 p += Xav[i] * (*si) * Xbv[i];
484 }
485 P(Xa.index2(),Xb.index2()) += p;
486 }
487 }
488}
489
490}//namespace detail
491
492
493/*
494 * prod_SPD - uBLAS Expression templates for X*X' and X*X'
495 * Functions are only defined for the type for which the operation is efficient
496 * ISSUE Although numerically symmetric, uBlas has no expression type to represent this property
497 * The result must be assigned to a symmetric container to exploit the symmetry
498 */
499
500template <class E1, class E2>
502{ // Provide ET result E1E2T_type of prod(matrix_expression<E1>,trans(matrix_expression<E2>)
503 typedef BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits<E2, ublas::scalar_identity<BOOST_UBLAS_TYPENAME E2::value_type> >::result_type E2T_type;
504 typedef BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits<BOOST_UBLAS_TYPENAME E1::value_type, E1,
505 BOOST_UBLAS_TYPENAME E2T_type::value_type, E2T_type>::result_type E1E2T_type;
506
507 // Provide ET result E1TE2_type of prod(trans(matrix_expression<E1>),matrix_expression<E2>)
508 typedef BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits<E1, ublas::scalar_identity<BOOST_UBLAS_TYPENAME E1::value_type> >::result_type E1T_type;
509 typedef BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits<BOOST_UBLAS_TYPENAME E1T_type::value_type, E1T_type,
510 BOOST_UBLAS_TYPENAME E2::value_type, E2>::result_type E1TE2_type;
511};
512
513
514template<class E> inline
516 prod_SPD (const ublas::matrix_expression<E>& X)
517/*
518 * Symmetric Positive (Semi) Definite product: X*X'
519 */
520{
521 // ISSUE: uBLAS post Boost 1_31_0 introduces a trans(const matrix_expression<E>& e) which propagates the const expression type
522 // Bypass this to avoid having to detect the Boost version
523 return prod( X, trans(const_cast<ublas::matrix_expression<E>&>(X) ));
524}
525
526template<class EX, class ES, class ET> inline
528 prod_SPD (const ublas::matrix_expression<EX>& X, const ublas::matrix_expression<ES>& S, ublas::matrix_expression<ET>& XStemp)
529/*
530 * Symmetric Positive (Semi) Definite product: X*(X*S)', XStemp = X*S
531 * Result symmetric if S is symmetric
532 */
533{
534 return prod( X, trans(prod(X,S,XStemp())) );
535}
536
537template<class EX, class ES, class ET> inline
538ET prod_SPD (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s, ublas::matrix_expression<ET>& Ptemp)
539/*
540 * Symmetric Positive (Semi) Definite product: X*diag_matrix(s)*X'
541 * Precond: Ptemp must be size conformant with the product
542 * DEPRECATED With explicit Ptemp, do not rely on it's value
543 */
544{
545 Vec::const_iterator si, send = s().end();
546 const EX& XX = X();
547 typename EX::const_iterator1 Xa = XX.begin1();
548 const typename EX::const_iterator1 Xend = XX.end1();
549 typename EX::const_iterator1 Xb;
550
551 // P(a,b) = sum(X.row(a) * s * X.row(b))
552 for (; Xa != Xend; ++Xa) // Iterate rows
553 {
554 typedef const ublas::matrix_row<const EX> EX_row;
555 EX_row Xav (ublas::row(XX, Xa.index1()));
556 Xb = Xa; // Start at the row Xa only one triangle of symmetric result required
557 for (; Xb != Xend; ++Xb)
558 {
559 typename EX::value_type p = 0; // Triple vector inner product
560 EX_row Xbv (ublas::row(XX, Xb.index1()));
561 for (si = s().begin(); si != send; ++si) {
562 Vec::size_type i = si.index();
563 p += Xav[i] * (*si) * Xbv[i];
564 }
565 Ptemp()(Xa.index1(),Xb.index1()) = p;
566 Ptemp()(Xb.index1(),Xa.index1()) = p;
567 }
568 }
569 return Ptemp();
570}
571
572template<class EX, class ES> inline
573SymMatrix prod_SPD (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s)
574/*
575 * Symmetric Positive (Semi) Definite product: X*diag_matrix(s)*X'
576 * TODO requires a prod_diag(X,s)
577 */
578{
579 SymMatrix Ptemp(X().size1(),X().size1());
580 Ptemp.clear();
581
582 (void) prod_SPD (X(),s(), Ptemp);
583 return Ptemp;
584}
585
586
587template<class E> inline
589 prod_SPDT (const ublas::matrix_expression<E>& X)
590/*
591 * Symmetric Positive (Semi) Definite product: X'*X
592 */
593{
594 // ISSUE: See prod_SPD
595 return prod( trans(const_cast<ublas::matrix_expression<E>&>(X) ), X);
596}
597
598template<class EX, class ES, class ET> inline
600 prod_SPDT (const ublas::matrix_expression<EX>& X, const ublas::matrix_expression<ES>& S, ublas::matrix_expression<ET>& SXtemp)
601/*
602 * Symmetric Positive (Semi) Definite product: (S*X)'*X, SXtemp = S*X
603 * Result symmetric if S is symmetric
604 */
605{
606 return prod( trans(prod(S,X,SXtemp())), X);
607}
608
609template<class EX, class ES, class ET> inline
610ET prod_SPDT (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s, ublas::matrix_expression<ET>& Ptemp)
611/*
612 * Symmetric Positive (Semi) Definite product: X'*diag_matrix(s)*X, Ptemp = return value
613 * Precond: Ptemp must be size conformant with the product
614 * DEPRECATED With explicit Ptemp, do not rely on it's value
615 */
616{
617 const EX& XX = X();
618 Vec::const_iterator si, send = s().end();
619 typename EX::const_iterator2 Xa = X.begin2();
620 const typename EX::const_iterator2 Xend = X.end2();
621 typename EX::const_iterator2 Xb;
622
623 // P(a,b) = sum(X.col(a) * s * X.col(b))
624 for (; Xa != Xend; ++Xa) // Iterate columns
625 {
626 typedef const ublas::matrix_column<const EX> EX_column;
627 EX_column Xav = ublas::column(XX, Xa.index2());
628 Xb = Xa; // Start at the column Xa only one triangle of symmetric result required
629 for (; Xb != Xend; ++Xb)
630 {
631 typename EX::value_type p = 0; // Triple vector inner product
632 EX_column Xbv = ublas::column(XX, Xb.index2());
633 for (si = s().begin(); si != send; ++si) {
634 Vec::size_type i = si.index();
635 p += Xav[i] * (*si) * Xbv[i];
636 }
637 Ptemp()(Xa.index2(),Xb.index2()) = p;
638 Ptemp()(Xb.index2(),Xa.index2()) = p;
639 }
640 }
641 return Ptemp();
642}
643
644template<class EX, class ES> inline
645SymMatrix prod_SPDT (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s)
646/*
647 * Symmetric Positive (Semi) Definite product: X'*diag_matrix(s)*X
648 * TODO requires a prod_diag(X,s)
649 */
650{
651 SymMatrix Ptemp(X().size1(),X().size1());
652 Ptemp.clear();
653
654 (void) prod_SPDT (X(),s(), Ptemp);
655 return Ptemp;
656}
657
658inline Vec::value_type prod_SPDT (const Vec& x, const Vec& s)
659/*
660 * Symmetric Positive (Semi) Definite product: x'*diag_matrix(s)*x
661 */
662{
663 const Vec::const_iterator xi_end = x.end();
664 Vec::const_iterator si = s.begin();
665 Vec::const_iterator xi=x.begin();
666 Vec::value_type p = Vec::value_type(0);
667 while (xi != xi_end)
668 {
669 p += (*xi) * (*si) * (*xi);
670 ++xi; ++ si;
671 }
672
673 return p;
674}
675
676inline Vec::value_type prod_SPDT (const Vec& x, const SymMatrix& S)
677/*
678 * Symmetric Positive (Semi) Definite multiply: p = x'*S*x
679 */
680{
681 return inner_prod (x, prod(S,x) );
682}
683
684}//namespace
BOOST_UBLAS_INLINE noalias_proxy(C &lval)
Definition uBLASmatrix.hpp:52
BOOST_UBLAS_INLINE void operator-=(const E &e)
Definition uBLASmatrix.hpp:69
BOOST_UBLAS_INLINE void operator=(const E &e)
Definition uBLASmatrix.hpp:57
BOOST_UBLAS_INLINE void operator+=(const E &e)
Definition uBLASmatrix.hpp:63
void mult_SPD(const MatrixX &X, const Vec &s, SymMatrix &P)
Definition uBLASmatrix.hpp:433
void mult_SPDT(const MatrixX &X, const Vec &s, SymMatrix &P)
Definition uBLASmatrix.hpp:462
Definition matSup.cpp:34
FMMatrix< detail::BaseDiagMatrix > DiagMatrix
Definition uBLASmatrix.hpp:330
FMMatrix< detail::BaseColMatrix > ColMatrix
Definition uBLASmatrix.hpp:326
FMMatrix< detail::BaseRowMatrix > RowMatrix
Definition uBLASmatrix.hpp:324
FMVec< detail::BaseDenseVector > DenseVec
Definition uBLASmatrix.hpp:333
FMMatrix< detail::SymMatrixWrapper< detail::BaseRowMatrix > > SymMatrix
Definition uBLASmatrix.hpp:327
FMVec< detail::BaseVector > Vec
Definition uBLASmatrix.hpp:323
prod_expression_result< E, E >::E1TE2_type prod_SPDT(const ublas::matrix_expression< E > &X)
Definition uBLASmatrix.hpp:589
FMMatrix< detail::SymMatrixWrapper< detail::BaseDenseRowMatrix > > DenseSymMatrix
Definition uBLASmatrix.hpp:337
EmptyTag
Definition uBLASmatrix.hpp:34
@ Empty
Definition uBLASmatrix.hpp:34
void identity(FMMatrix< Base > &I)
Definition uBLASmatrix.hpp:412
FMMatrix< detail::BaseDenseLowerTriMatrix > DenseLTriMatrix
Definition uBLASmatrix.hpp:339
ublas::matrix_vector_range< FMMatrix< Base > > diag(FMMatrix< Base > &M, std::size_t n)
Definition uBLASmatrix.hpp:383
BOOST_UBLAS_INLINE detail::noalias_proxy< E > noalias(ublas::matrix_expression< E > &lvalue)
Definition uBLASmatrix.hpp:85
const ublas::triangular_adaptor< const M, ublas::upper > UpperTri(const M &m)
Definition uBLASmatrix.hpp:357
FMMatrix< detail::BaseDenseUpperTriMatrix > DenseUTriMatrix
Definition uBLASmatrix.hpp:338
FMMatrix< detail::BaseUpperTriMatrix > UTriMatrix
Definition uBLASmatrix.hpp:328
RowMatrix Matrix
Definition uBLASmatrix.hpp:325
FMMatrix< detail::BaseLowerTriMatrix > LTriMatrix
Definition uBLASmatrix.hpp:329
FMMatrix< detail::BaseDenseDiagMatrix > DenseDiagMatrix
Definition uBLASmatrix.hpp:340
FMMatrix< detail::BaseDenseRowMatrix > DenseRowMatrix
Definition uBLASmatrix.hpp:334
const ublas::triangular_adaptor< const M, ublas::lower > LowerTri(const M &m)
Definition uBLASmatrix.hpp:368
FMMatrix< detail::BaseDenseColMatrix > DenseColMatrix
Definition uBLASmatrix.hpp:336
prod_expression_result< E, E >::E1E2T_type prod_SPD(const ublas::matrix_expression< E > &X)
Definition uBLASmatrix.hpp:516
DenseRowMatrix DenseMatrix
Definition uBLASmatrix.hpp:335
BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits< BOOST_UBLAS_TYPENAMEE1T_type::value_type, E1T_type, BOOST_UBLAS_TYPENAMEE2::value_type, E2 >::result_type E1TE2_type
Definition uBLASmatrix.hpp:510
BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits< BOOST_UBLAS_TYPENAMEE1::value_type, E1, BOOST_UBLAS_TYPENAMEE2T_type::value_type, E2T_type >::result_type E1E2T_type
Definition uBLASmatrix.hpp:505
BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits< E1, ublas::scalar_identity< BOOST_UBLAS_TYPENAMEE1::value_type > >::result_type E1T_type
Definition uBLASmatrix.hpp:508
BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits< E2, ublas::scalar_identity< BOOST_UBLAS_TYPENAMEE2::value_type > >::result_type E2T_type
Definition uBLASmatrix.hpp:503
#define BOOST_UBLAS_TYPENAME
Definition uBLASmatrix.hpp:38