invert_psd_matrix.h 3.3 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677
  1. // Ceres Solver - A fast non-linear least squares minimizer
  2. // Copyright 2023 Google Inc. All rights reserved.
  3. // http://ceres-solver.org/
  4. //
  5. // Redistribution and use in source and binary forms, with or without
  6. // modification, are permitted provided that the following conditions are met:
  7. //
  8. // * Redistributions of source code must retain the above copyright notice,
  9. // this list of conditions and the following disclaimer.
  10. // * Redistributions in binary form must reproduce the above copyright notice,
  11. // this list of conditions and the following disclaimer in the documentation
  12. // and/or other materials provided with the distribution.
  13. // * Neither the name of Google Inc. nor the names of its contributors may be
  14. // used to endorse or promote products derived from this software without
  15. // specific prior written permission.
  16. //
  17. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18. // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19. // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20. // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21. // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22. // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23. // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24. // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25. // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26. // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27. // POSSIBILITY OF SUCH DAMAGE.
  28. //
  29. // Author: sameeragarwal@google.com (Sameer Agarwal)
  30. #ifndef CERES_INTERNAL_INVERT_PSD_MATRIX_H_
  31. #define CERES_INTERNAL_INVERT_PSD_MATRIX_H_
  32. #include "Eigen/Dense"
  33. #include "ceres/internal/eigen.h"
  34. #include "glog/logging.h"
  35. namespace ceres::internal {
  36. // Helper routine to compute the inverse or pseudo-inverse of a
  37. // symmetric positive semi-definite matrix.
  38. //
  39. // assume_full_rank controls whether a Cholesky factorization or an
  40. // Singular Value Decomposition is used to compute the inverse and the
  41. // pseudo-inverse respectively.
  42. //
  43. // The template parameter kSize can either be Eigen::Dynamic or a
  44. // positive integer equal to the number of rows of m.
  45. template <int kSize>
  46. typename EigenTypes<kSize, kSize>::Matrix InvertPSDMatrix(
  47. const bool assume_full_rank,
  48. const typename EigenTypes<kSize, kSize>::Matrix& m) {
  49. using MType = typename EigenTypes<kSize, kSize>::Matrix;
  50. const int size = m.rows();
  51. // If the matrix can be assumed to be full rank, then if it is small
  52. // (< 5) and fixed size, use Eigen's optimized inverse()
  53. // implementation.
  54. //
  55. // https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html#title3
  56. if (assume_full_rank) {
  57. if (kSize > 0 && kSize < 5) {
  58. return m.inverse();
  59. }
  60. return m.template selfadjointView<Eigen::Upper>().llt().solve(
  61. MType::Identity(size, size));
  62. }
  63. // For a thin SVD the number of columns of the matrix need to be dynamic.
  64. using SVDMType = typename EigenTypes<kSize, Eigen::Dynamic>::Matrix;
  65. Eigen::JacobiSVD<SVDMType> svd(m, Eigen::ComputeThinU | Eigen::ComputeThinV);
  66. return svd.solve(MType::Identity(size, size));
  67. }
  68. } // namespace ceres::internal
  69. #endif // CERES_INTERNAL_INVERT_PSD_MATRIX_H_