Belle II Software development
matrix.test.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9/*
10This file contains some basic investigation on numerically matrix algorithms.
11*/
12
13#include <gtest/gtest.h>
14
15#include <vector>
16#include <Eigen/Dense>
17
18
19// 1. Study of matrix inversions
20
21// TL;DR use a rank revealing matrix inversion (usually involving pivoting)
22// * pivoted householder QR decomposition
23// * full pivoted LU
24
25// no good for inversion if rank of the matrix is not full.
26// TEST(TrackFindingCDCTest, TMatrixDSym_invert)
27// {
28// TMatrixDSym weight(2);
29// weight(0,0) = 2;
30// weight(0,1) = 0;
31// weight(1,0) = 0;
32// weight(1,1) = 0;
33
34// EXPECT_EQ(2, weight.GetNrows());
35// EXPECT_EQ(2, weight.GetNcols());
36
37// TMatrixDSym cov = weight;
38// cov.Invert();
39
40// EXPECT_EQ(2, weight.GetNrows());
41// EXPECT_EQ(2, weight.GetNcols());
42
43// // Already fails here
44// EXPECT_EQ(2, cov.GetNrows());
45// EXPECT_EQ(2, cov.GetNcols());
46
47// EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
48// EXPECT_NEAR(0.0, cov(0,1), 10e-7);
49// EXPECT_NEAR(0.0, cov(1,0), 10e-7);
50// EXPECT_NEAR(0.0, cov(1,1), 10e-7);
51// }
52
53// no good for inversion if rank of the matrix is not full.
54// TEST(TrackFindingCDCTest, SMatrixDSym_inverse)
55// {
56// using CovarianceMatrix = ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2>
57// >;
58// using WeightMatrix = ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2> >;
59
60// WeightMatrix weight;
61// weight(0,0) = 2;
62// weight(0,1) = 0;
63// weight(1,0) = 0;
64// weight(1,1) = 0;
65
66// int state = 0;
67// CovarianceMatrix cov = weight.Inverse(state);
68
69// EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
70// EXPECT_NEAR(0.0, cov(0,1), 10e-7);
71// EXPECT_NEAR(0.0, cov(1,0), 10e-7);
72// EXPECT_NEAR(0.0, cov(1,1), 10e-7);
73// }
74
75// no good for inversion if rank of the matrix is not full.
76// TEST(TrackFindingCDCTest, SMatrixDSym_inverseChol)
77// {
78// using CovarianceMatrix = ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2>
79// >;
80// using WeightMatrix = ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2> >;
81
82// WeightMatrix weight;
83// weight(0,0) = 2;
84// weight(0,1) = 0;
85// weight(1,0) = 0;
86// weight(1,1) = 0;
87
88// int state = 0;
89// CovarianceMatrix cov = weight.InverseChol(state);
90
91// EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
92// EXPECT_NEAR(0.0, cov(0,1), 10e-7);
93// EXPECT_NEAR(0.0, cov(1,0), 10e-7);
94// EXPECT_NEAR(0.0, cov(1,1), 10e-7);
95// }
96
97// no good for inversion if rank of the matrix is not full.
98// TEST(TrackFindingCDCTest, Eigen_inverse)
99// {
100// Eigen::Matrix<double, 2, 2> weight;
101// weight(0,0) = 2;
102// weight(0,1) = 0;
103// weight(1,0) = 0;
104// weight(1,1) = 0;
105
106// Eigen::Matrix<double, 2, 2> cov = weight.inverse();
107
108// EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
109// EXPECT_NEAR(0.0, cov(0,1), 10e-7);
110// EXPECT_NEAR(0.0, cov(1,0), 10e-7);
111// EXPECT_NEAR(0.0, cov(1,1), 10e-7);
112// }
113
114// no good for inversion if rank of the matrix is not full.
115// TEST(TrackFindingCDCTest, Eigen_llt_inverse)
116// {
117// Eigen::Matrix<double, 2, 2> weight;
118// weight(0,0) = 2;
119// weight(0,1) = 0;
120// weight(1,0) = 0;
121// weight(1,1) = 0;
122
123// Eigen::Matrix<double, 2, 2> cov = weight.llt().inverse();
124
125// EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
126// EXPECT_NEAR(0.0, cov(0,1), 10e-7);
127// EXPECT_NEAR(0.0, cov(1,0), 10e-7);
128// EXPECT_NEAR(0.0, cov(1,1), 10e-7);
129// }
130
131// No good for inversion if rank of the matrix is not full.
132// TEST(TrackFindingCDCTest, Eigen_lldt_inverse)
133// {
134// Eigen::Matrix<double, 2, 2> weight;
135// weight(0,0) = 2;
136// weight(0,1) = 0;
137// weight(1,0) = 0;
138// weight(1,1) = 0;
139
140// Eigen::Matrix<double, 2, 2> cov = weight.lldt().inverse();
141
142// EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
143// EXPECT_NEAR(0.0, cov(0,1), 10e-7);
144// EXPECT_NEAR(0.0, cov(1,0), 10e-7);
145// EXPECT_NEAR(0.0, cov(1,1), 10e-7);
146// }
147
148TEST(TrackFindingCDCTest, Eigen_fullPivLu_inverse)
149{
150 Eigen::Matrix<double, 2, 2> weight;
151 weight(0, 0) = 2;
152 weight(0, 1) = 0;
153 weight(1, 0) = 0;
154 weight(1, 1) = 0;
155
156 Eigen::Matrix<double, 2, 2> cov = weight.fullPivLu().inverse();
157
158 EXPECT_NEAR(1.0 / 2.0, cov(0, 0), 10e-7);
159 EXPECT_NEAR(0.0, cov(0, 1), 10e-7);
160 EXPECT_NEAR(0.0, cov(1, 0), 10e-7);
161 EXPECT_NEAR(0.0, cov(1, 1), 10e-7);
162}
163
164// No good for inversion if rank of the matrix is not full.
165// TEST(TrackFindingCDCTest, Eigen_partialPivLu_inverse)
166// {
167// Eigen::Matrix<double, 2, 2> weight;
168// weight(0,0) = 2;
169// weight(0,1) = 0;
170// weight(1,0) = 0;
171// weight(1,1) = 0;
172
173// Eigen::Matrix<double, 2, 2> cov = weight.partialPivLu().inverse();
174
175// EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
176// EXPECT_NEAR(0.0, cov(0,1), 10e-7);
177// EXPECT_NEAR(0.0, cov(1,0), 10e-7);
178// EXPECT_NEAR(0.0, cov(1,1), 10e-7);
179// }
180
181// Can not do inversion
182// TEST(TrackFindingCDCTest, Eigen_householderQr_inverse)
183// {
184// Eigen::Matrix<double, 2, 2> weight;
185// weight(0,0) = 2;
186// weight(0,1) = 0;
187// weight(1,0) = 0;
188// weight(1,1) = 0;
189
190// Eigen::Matrix<double, 2, 2> cov = weight.householderQr().inverse();
191
192// EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
193// EXPECT_NEAR(0.0, cov(0,1), 10e-7);
194// EXPECT_NEAR(0.0, cov(1,0), 10e-7);
195// EXPECT_NEAR(0.0, cov(1,1), 10e-7);
196// }
197
198/* Method of choice for covariance <-> weight matrices */
199TEST(TrackFindingCDCTest, Eigen_colPivHouseholderQr_inverse)
200{
201 Eigen::Matrix<double, 2, 2> weight;
202 weight(0, 0) = 2;
203 weight(0, 1) = 0;
204 weight(1, 0) = 0;
205 weight(1, 1) = 0;
206
207 Eigen::Matrix<double, 2, 2> cov = weight.colPivHouseholderQr().inverse();
208
209 EXPECT_NEAR(1.0 / 2.0, cov(0, 0), 10e-7);
210 EXPECT_NEAR(0.0, cov(0, 1), 10e-7);
211 EXPECT_NEAR(0.0, cov(1, 0), 10e-7);
212 EXPECT_NEAR(0.0, cov(1, 1), 10e-7);
213}
214
215TEST(TrackFindingCDCTest, Eigen_fullPivHouseholderQr_inverse)
216{
217 Eigen::Matrix<double, 2, 2> weight;
218 weight(0, 0) = 2;
219 weight(0, 1) = 0;
220 weight(1, 0) = 0;
221 weight(1, 1) = 0;
222
223 Eigen::Matrix<double, 2, 2> cov = weight.fullPivHouseholderQr().inverse();
224
225 EXPECT_NEAR(1.0 / 2.0, cov(0, 0), 10e-7);
226 EXPECT_NEAR(0.0, cov(0, 1), 10e-7);
227 EXPECT_NEAR(0.0, cov(1, 0), 10e-7);
228 EXPECT_NEAR(0.0, cov(1, 1), 10e-7);
229}