Added the inv_upper_triangular() and inv_upper_triangular()

functions.

--HG--
extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%402402

dlib/matrix/matrix_utilities.h

@@ -2615,6 +2615,93 @@ namespace dlib
         const matrix_exp<EXP>& m
     ) { return inv_helper<EXP,matrix_exp<EXP>::NR>::inv(m); }
 
+// ----------------------------------------------------------------------------------------
+
+    template <typename EXP>
+    const typename matrix_exp<EXP>::matrix_type  inv_lower_triangular (
+        const matrix_exp<EXP>& A 
+    )
+    {
+        DLIB_ASSERT(A.nr() == A.nc(), 
+            "\tconst matrix inv_lower_triangular(const matrix_exp& A)"
+            << "\n\tA must be a square matrix"
+            << "\n\tA.nr(): " << A.nr()
+            << "\n\tA.nc(): " << A.nc() 
+            );
+
+        typedef typename matrix_exp<EXP>::matrix_type matrix_type;
+        typedef typename matrix_type::type type;
+
+        matrix_type m(A);
+
+        for(long c = 0; c < m.nc(); ++c)
+        {
+            if( m(c,c) == 0 )
+            {
+                // there isn't an inverse so just give up
+                return m;
+            }
+
+            // compute m(c,c)
+            m(c,c) = 1/m(c,c);
+
+            // compute the values in column c that are below m(c,c).
+            // We do this by just doing the same thing we do for upper triangular
+            // matrices because we take the transpose of m which turns m into an
+            // upper triangular matrix.
+            for(long r = 0; r < c; ++r)
+            {
+                const long n = c-r;
+                m(c,r) = -m(c,c)*subm(trans(m),r,r,1,n)*subm(trans(m),r,c,n,1);
+            }
+        }
+
+        return m;
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename EXP>
+    const typename matrix_exp<EXP>::matrix_type  inv_upper_triangular (
+        const matrix_exp<EXP>& A 
+    )
+    {
+        DLIB_ASSERT(A.nr() == A.nc(), 
+            "\tconst matrix inv_upper_triangular(const matrix_exp& A)"
+            << "\n\tA must be a square matrix"
+            << "\n\tA.nr(): " << A.nr()
+            << "\n\tA.nc(): " << A.nc() 
+            );
+
+        typedef typename matrix_exp<EXP>::matrix_type matrix_type;
+        typedef typename matrix_type::type type;
+
+        matrix_type m(A);
+
+        for(long c = 0; c < m.nc(); ++c)
+        {
+            if( m(c,c) == 0 )
+            {
+                // there isn't an inverse so just give up
+                return m;
+            }
+
+            // compute m(c,c)
+            m(c,c) = 1/m(c,c);
+
+            // compute the values in column c that are above m(c,c)
+            for(long r = 0; r < c; ++r)
+            {
+                const long n = c-r;
+                m(r,c) = -m(c,c)*subm(m,r,r,1,n)*subm(m,r,c,n,1);
+            }
+        }
+
+        return m;
+
+    }
+
 // ----------------------------------------------------------------------------------------
 
     template <

dlib/matrix/matrix_utilities_abstract.h

@@ -679,6 +679,38 @@ namespace dlib
                   will have a bogus value.  I.e. it won't be a decomposition.
     !*/
 
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::matrix_type inv_lower_triangular (
+        const matrix_exp& A
+    );
+    /*!
+        requires
+            - A is a square matrix
+        ensures
+            - if (A is lower triangular) then
+                - returns the inverse of A. 
+            - else
+                - returns a matrix with the same dimensions as A but it 
+                  will have a bogus value.  I.e. it won't be an inverse.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::matrix_type inv_upper_triangular (
+        const matrix_exp& A
+    );
+    /*!
+        requires
+            - A is a square matrix
+        ensures
+            - if (A is upper triangular) then
+                - returns the inverse of A. 
+            - else
+                - returns a matrix with the same dimensions as A but it 
+                  will have a bogus value.  I.e. it won't be an inverse.
+    !*/
+
 // ----------------------------------------------------------------------------------------
 // ----------------------------------------------------------------------------------------
 //                              Statistics

dlib/test/matrix.cpp

@@ -1302,6 +1302,82 @@ namespace
 
         }
 
+
+        {
+
+            const double stuff[] = { 
+                1, 2, 3,
+                6, 3, 3, 
+                7, 3, 9};
+
+            matrix<double,3,3> m(stuff);
+
+            // make m be symmetric
+            m = m*trans(m);
+
+            matrix<double> L = cholesky_decomposition(m);
+            DLIB_CASSERT(equal(L*trans(L), m), "");
+
+            DLIB_CASSERT(equal(inv(m), inv_upper_triangular(trans(L))*inv_lower_triangular((L))), "") 
+            DLIB_CASSERT(equal(round_zeros(inv_upper_triangular(trans(L))*trans(L),1e-10), identity_matrix<double>(3), 1e-10),""); 
+            DLIB_CASSERT(equal(round_zeros(inv_lower_triangular((L))*(L),1e-10) ,identity_matrix<double>(3),1e-10),""); 
+
+        }
+
+        {
+
+            const double stuff[] = { 
+                1, 2, 3, 6, 3, 4,
+                6, 3, 3, 1, 2, 3,
+                7, 3, 9, 54.3, 5, 3,
+                -6, 3, -3, 1, 2, 3,
+                1, 2, 3, 5, -3, 3,
+                7, 3, -9, 54.3, 5, 3
+                };
+
+            matrix<double,6,6> m(stuff);
+
+            // make m be symmetric
+            m = m*trans(m);
+
+            matrix<double> L = cholesky_decomposition(m);
+            DLIB_CASSERT(equal(L*trans(L), m, 1e-10), L*trans(L)-m);
+
+            DLIB_CASSERT(equal(inv(m), inv_upper_triangular(trans(L))*inv_lower_triangular((L))), "") 
+            DLIB_CASSERT(equal(inv(m), trans(inv_lower_triangular(L))*inv_lower_triangular((L))), "") 
+            DLIB_CASSERT(equal(inv(m), trans(inv_lower_triangular(L))*trans(inv_upper_triangular(trans(L)))), "") 
+            DLIB_CASSERT(equal(round_zeros(inv_upper_triangular(trans(L))*trans(L),1e-10) , identity_matrix<double>(6), 1e-10),
+                         round_zeros(inv_upper_triangular(trans(L))*trans(L),1e-10)); 
+            DLIB_CASSERT(equal(round_zeros(inv_lower_triangular((L))*(L),1e-10) ,identity_matrix<double>(6), 1e-10),
+                         round_zeros(inv_lower_triangular((L))*(L),1e-10)); 
+
+        }
+
+        {
+
+            matrix<double,6,6> m(identity_matrix<double>(6)*4.5);
+
+            matrix<double> L = cholesky_decomposition(m);
+            DLIB_CASSERT(equal(L*trans(L), m, 1e-10), L*trans(L)-m);
+
+            DLIB_CASSERT(equal(inv(m), inv_upper_triangular(trans(L))*inv_lower_triangular((L))), "") 
+            DLIB_CASSERT(equal(round_zeros(inv_upper_triangular(trans(L))*trans(L),1e-10) , identity_matrix<double>(6), 1e-10),
+                         round_zeros(inv_upper_triangular(trans(L))*trans(L),1e-10)); 
+            DLIB_CASSERT(equal(round_zeros(inv_lower_triangular((L))*(L),1e-10) ,identity_matrix<double>(6), 1e-10),
+                         round_zeros(inv_lower_triangular((L))*(L),1e-10)); 
+
+        }
+
+        {
+
+            matrix<double,6,6> m(identity_matrix<double>(6)*4.5);
+            m(1,4) = 2;
+
+            DLIB_CASSERT(dlib::equal(inv_upper_triangular(m), inv(m),1e-10), inv_upper_triangular(m)-inv(m));
+            DLIB_CASSERT(dlib::equal(inv_lower_triangular(trans(m)), inv(trans(m)),1e-10), inv_lower_triangular(trans(m))-inv(trans(m)));
+
+        }
+
     }