1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
|
// Copyright (C) 2009 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
// This code was adapted from code from the JAMA part of NIST's TNT library.
// See: http://math.nist.gov/tnt/
#ifndef DLIB_MATRIX_CHOLESKY_DECOMPOSITION_H
#define DLIB_MATRIX_CHOLESKY_DECOMPOSITION_H
#include "matrix.h"
#include "matrix_utilities.h"
#include "matrix_subexp.h"
#include <cmath>
#ifdef DLIB_USE_LAPACK
#include "lapack/potrf.h"
#endif
#include "matrix_trsm.h"
namespace dlib
{
template <
typename matrix_exp_type
>
class cholesky_decomposition
{
public:
const static long NR = matrix_exp_type::NR;
const static long NC = matrix_exp_type::NC;
typedef typename matrix_exp_type::type type;
typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
typedef typename matrix_exp_type::layout_type layout_type;
typedef matrix<type,0,0,mem_manager_type,layout_type> matrix_type;
typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;
// You have supplied an invalid type of matrix_exp_type. You have
// to use this object with matrices that contain float or double type data.
COMPILE_TIME_ASSERT((is_same_type<float, type>::value ||
is_same_type<double, type>::value ));
template <typename EXP>
cholesky_decomposition(
const matrix_exp<EXP>& A
);
bool is_spd(
) const;
const matrix_type& get_l(
) const;
template <typename EXP>
const typename EXP::matrix_type solve (
const matrix_exp<EXP>& B
) const;
private:
matrix_type L_; // lower triangular factor
bool isspd; // true if matrix to be factored was SPD
};
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// Member functions
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
template <typename matrix_exp_type>
bool cholesky_decomposition<matrix_exp_type>::
is_spd(
) const
{
return isspd;
}
// ----------------------------------------------------------------------------------------
template <typename matrix_exp_type>
const typename cholesky_decomposition<matrix_exp_type>::matrix_type& cholesky_decomposition<matrix_exp_type>::
get_l(
) const
{
return L_;
}
// ----------------------------------------------------------------------------------------
template <typename matrix_exp_type>
template <typename EXP>
cholesky_decomposition<matrix_exp_type>::
cholesky_decomposition(
const matrix_exp<EXP>& A_
)
{
using std::sqrt;
COMPILE_TIME_ASSERT((is_same_type<type, typename EXP::type>::value));
// make sure requires clause is not broken
DLIB_ASSERT(A_.nr() == A_.nc() && A_.size() > 0,
"\tcholesky_decomposition::cholesky_decomposition(A_)"
<< "\n\tYou can only use this on square matrices"
<< "\n\tA_.nr(): " << A_.nr()
<< "\n\tA_.nc(): " << A_.nc()
<< "\n\tA_.size(): " << A_.size()
<< "\n\tthis: " << this
);
#ifdef DLIB_USE_LAPACK
L_ = A_;
const type eps = max(abs(diag(L_)))*std::sqrt(std::numeric_limits<type>::epsilon())/100;
// check if the matrix is actually symmetric
bool is_symmetric = true;
for (long r = 0; r < L_.nr() && is_symmetric; ++r)
{
for (long c = r+1; c < L_.nc() && is_symmetric; ++c)
{
// this is approximately doing: is_symmetric = is_symmetric && ( L_(k,j) == L_(j,k))
is_symmetric = is_symmetric && (std::abs(L_(r,c) - L_(c,r)) < eps );
}
}
// now compute the actual cholesky decomposition
int info = lapack::potrf('L', L_);
// check if it's really SPD
if (info == 0 && is_symmetric && min(abs(diag(L_))) > eps*100)
isspd = true;
else
isspd = false;
L_ = lowerm(L_);
#else
const_temp_matrix<EXP> A(A_);
isspd = true;
const long n = A.nc();
L_.set_size(n,n);
const type eps = max(abs(diag(A)))*std::sqrt(std::numeric_limits<type>::epsilon())/100;
// Main loop.
for (long j = 0; j < n; j++)
{
type d(0.0);
for (long k = 0; k < j; k++)
{
type s(0.0);
for (long i = 0; i < k; i++)
{
s += L_(k,i)*L_(j,i);
}
// if L_(k,k) != 0
if (std::abs(L_(k,k)) > eps)
{
s = (A(j,k) - s)/L_(k,k);
}
else
{
s = (A(j,k) - s);
isspd = false;
}
L_(j,k) = s;
d = d + s*s;
// this is approximately doing: isspd = isspd && ( A(k,j) == A(j,k))
isspd = isspd && (std::abs(A(k,j) - A(j,k)) < eps );
}
d = A(j,j) - d;
isspd = isspd && (d > eps);
L_(j,j) = sqrt(d > 0.0 ? d : 0.0);
for (long k = j+1; k < n; k++)
{
L_(j,k) = 0.0;
}
}
#endif
}
// ----------------------------------------------------------------------------------------
template <typename matrix_exp_type>
template <typename EXP>
const typename EXP::matrix_type cholesky_decomposition<matrix_exp_type>::
solve(
const matrix_exp<EXP>& B
) const
{
COMPILE_TIME_ASSERT((is_same_type<type, typename EXP::type>::value));
// make sure requires clause is not broken
DLIB_ASSERT(L_.nr() == B.nr(),
"\tconst matrix cholesky_decomposition::solve(B)"
<< "\n\tInvalid arguments were given to this function."
<< "\n\tL_.nr(): " << L_.nr()
<< "\n\tB.nr(): " << B.nr()
<< "\n\tthis: " << this
);
matrix<type, NR, EXP::NC, mem_manager_type, layout_type> X(B);
using namespace blas_bindings;
// Solve L*y = b;
triangular_solver(CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, L_, X);
// Solve L'*X = Y;
triangular_solver(CblasLeft, CblasLower, CblasTrans, CblasNonUnit, L_, X);
return X;
}
// ----------------------------------------------------------------------------------------
}
#endif // DLIB_MATRIX_CHOLESKY_DECOMPOSITION_H
|
