ACloudViewer  3.9.4
A Modern Library for 3D Data Processing
least_absolute_deviations_test.cc
Go to the documentation of this file.
1 // Copyright (c) 2018, ETH Zurich and UNC Chapel Hill.
2 // All rights reserved.
3 //
4 // Redistribution and use in source and binary forms, with or without
5 // modification, are permitted provided that the following conditions are met:
6 //
7 // * Redistributions of source code must retain the above copyright
8 // notice, this list of conditions and the following disclaimer.
9 //
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 //
14 // * Neither the name of ETH Zurich and UNC Chapel Hill nor the names of
15 // its contributors may be used to endorse or promote products derived
16 // from this software without specific prior written permission.
17 //
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
22 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28 // POSSIBILITY OF SUCH DAMAGE.
29 //
30 // Author: Johannes L. Schoenberger (jsch-at-demuc-dot-de)
31 
32 #define TEST_NAME "optim/least_absolute_deviations"
33 #include "util/testing.h"
34 
35 #include <Eigen/Dense>
36 
38 #include "util/random.h"
39 
40 using namespace colmap;
41 
42 BOOST_AUTO_TEST_CASE(TestOverDetermined) {
43  Eigen::SparseMatrix<double> A(4, 3);
44  for (int i = 0; i < A.rows(); ++i) {
45  for (int j = 0; j < A.cols(); ++j) {
46  A.insert(i, j) = i * A.cols() + j + 1;
47  }
48  }
49  A.coeffRef(0, 0) = 10;
50 
51  Eigen::VectorXd b(A.rows());
52  for (int i = 0; i < b.size(); ++i) {
53  b(i) = i + 1;
54  }
55 
56  Eigen::VectorXd x = Eigen::VectorXd::Zero(A.cols());
57 
59  BOOST_CHECK(SolveLeastAbsoluteDeviations(options, A, b, &x));
60 
61  // Reference solution obtained with Boyd's Matlab implementation.
62  const Eigen::Vector3d x_ref(0, 0, 1 / 3.0);
63  BOOST_CHECK(x.isApprox(x_ref));
64 
65  const Eigen::VectorXd residual = A * x - b;
66  BOOST_CHECK_LE(residual.norm(), 1e-6);
67 }
68 
69 BOOST_AUTO_TEST_CASE(TestWellDetermined) {
70  Eigen::SparseMatrix<double> A(3, 3);
71  for (int i = 0; i < A.rows(); ++i) {
72  for (int j = 0; j < A.cols(); ++j) {
73  A.insert(i, j) = i * A.cols() + j + 1;
74  }
75  }
76  A.coeffRef(0, 0) = 10;
77 
78  Eigen::VectorXd b(A.rows());
79  for (int i = 0; i < b.size(); ++i) {
80  b(i) = i + 1;
81  }
82 
83  Eigen::VectorXd x = Eigen::VectorXd::Zero(A.cols());
84 
86  BOOST_CHECK(SolveLeastAbsoluteDeviations(options, A, b, &x));
87 
88  // Reference solution obtained with Boyd's Matlab implementation.
89  const Eigen::Vector3d x_ref(0, 0, 1 / 3.0);
90  BOOST_CHECK(x.isApprox(x_ref));
91 
92  const Eigen::VectorXd residual = A * x - b;
93  BOOST_CHECK_LE(residual.norm(), 1e-6);
94 }
95 
96 BOOST_AUTO_TEST_CASE(TestUnderDetermined) {
97  // In this case, the system is rank-deficient and not positive semi-definite.
98  Eigen::SparseMatrix<double> A(2, 3);
99  Eigen::VectorXd b(A.rows());
100  Eigen::VectorXd x = Eigen::VectorXd::Zero(A.cols());
102  BOOST_CHECK(!SolveLeastAbsoluteDeviations(options, A, b, &x));
103 }
const double * e
BOOST_AUTO_TEST_CASE(TestOverDetermined)
normal_z x
bool SolveLeastAbsoluteDeviations(const LeastAbsoluteDeviationsOptions &options, const Eigen::SparseMatrix< double > &A, const Eigen::VectorXd &b, Eigen::VectorXd *x)