Added QR, LU, cholesky, and eigenvalue decomposition classes.

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

dlib/matrix.h

@@ -8,6 +8,7 @@
 #include "matrix/matrix_subexp.h"
 #include "matrix/matrix_math_functions.h"
 #include "matrix/matrix_assign.h"
+#include "matrix/matrix_la.h"
 
 #ifdef DLIB_USE_BLAS
 #include "matrix/matrix_blas_bindings.h"

dlib/matrix/matrix_cholesky.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.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>
+
+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 typename matrix_exp_type::matrix_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 matrix<type,matrix_exp_type::NR,EXP::NC,mem_manager_type,layout_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
+            );
+
+        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;
+            }
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    template <typename EXP>
+    const matrix<typename matrix_exp_type::type,
+                 matrix_exp_type::NR, 
+                 EXP::NC, 
+                 typename matrix_exp_type::mem_manager_type, 
+                 typename matrix_exp_type::layout_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); 
+
+        const long nx = B.nc();
+
+        const long n = L_.nr();
+
+        // Solve L*y = b;
+        for (long j=0; j< nx; j++)
+        {
+            for (long k = 0; k < n; k++) 
+            {
+                for (long i = 0; i < k; i++) 
+                    X(k,j) -= X(i,j)*L_(k,i);
+                X(k,j) /= L_(k,k);
+            }
+        }
+
+        // Solve L'*X = Y;
+        for (long j=0; j<nx; j++)
+        {
+            for (long k = n-1; k >= 0; k--) 
+            {
+                for (long i = k+1; i < n; i++) 
+                    X(k,j) -= X(i,j)*L_(i,k);
+                X(k,j) /= L_(k,k);
+            }
+        }
+
+        return X;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+
+
+} 
+
+#endif // DLIB_MATRIX_CHOLESKY_DECOMPOSITION_H 
+
+
+
+

dlib/matrix/matrix_eigenvalue.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.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_EIGENVALUE_DECOMPOSITION_H
+#define DLIB_MATRIX_EIGENVALUE_DECOMPOSITION_H
+
+#include "matrix.h" 
+#include "matrix_utilities.h"
+#include "matrix_subexp.h"
+#include <algorithm>
+#include <complex>
+#include <cmath>
+
+namespace dlib 
+{
+
+    template <
+        typename matrix_exp_type
+        >
+    class eigenvalue_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 typename matrix_exp_type::matrix_type matrix_type;
+        typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;
+
+        typedef matrix<std::complex<type>,0,0,mem_manager_type,layout_type> complex_matrix_type;
+        typedef matrix<std::complex<type>,NR,1,mem_manager_type,layout_type> complex_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>
+        eigenvalue_decomposition(
+            const matrix_exp<EXP>& A
+        ); 
+
+        long dim (
+        ) const;
+
+        const complex_column_vector_type get_eigenvalues (
+        ) const;
+
+        const column_vector_type& get_real_eigenvalues (
+        ) const;
+
+        const column_vector_type& get_imag_eigenvalues (
+        ) const;
+
+        const complex_matrix_type get_v (
+        ) const;
+
+        const complex_matrix_type get_d (
+        ) const; 
+
+        const matrix_type& get_pseudo_v (
+        ) const;
+
+        const matrix_type get_pseudo_d (
+        ) const; 
+
+    private:
+
+        /** Row and column dimension (square matrix).  */
+        int n;
+
+        bool issymmetric; 
+
+        /** Arrays for internal storage of eigenvalues. */
+
+        column_vector_type d;         /* real part */
+        column_vector_type e;         /* img part */
+
+        /** Array for internal storage of eigenvectors. */
+        matrix_type V;
+
+        /** Array for internal storage of nonsymmetric Hessenberg form.
+        @serial internal storage of nonsymmetric Hessenberg form.
+        */
+        matrix_type H;
+
+
+        /** Working storage for nonsymmetric algorithm.
+        @serial working storage for nonsymmetric algorithm.
+        */
+        column_vector_type ort;
+
+        // Symmetric Householder reduction to tridiagonal form.
+        void tred2();
+
+
+        // Symmetric tridiagonal QL algorithm.
+        void tql2 ();
+
+
+        // Nonsymmetric reduction to Hessenberg form.
+        void orthes ();
+
+
+        // Complex scalar division.
+        type cdivr, cdivi;
+        void cdiv(type xr, type xi, type yr, type yi);
+
+
+        // Nonsymmetric reduction from Hessenberg to real Schur form.
+        void hqr2 (); 
+    };
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+//                        Public member functions
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    template <typename EXP>
+    eigenvalue_decomposition<matrix_exp_type>::
+    eigenvalue_decomposition(
+        const matrix_exp<EXP>& A__
+    ) 
+    {
+        COMPILE_TIME_ASSERT((is_same_type<type, typename EXP::type>::value));
+
+
+        const_temp_matrix<EXP> A(A__);
+
+        // make sure requires clause is not broken
+        DLIB_ASSERT(A.nr() == A.nc() && A.size() > 0,
+            "\teigenvalue_decomposition::eigenvalue_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
+            );
+
+
+        n = A.nc();
+        V.set_size(n,n);
+        d.set_size(n);
+        e.set_size(n);
+
+
+        issymmetric = true;
+        for (int j = 0; (j < n) && issymmetric; j++) 
+        {
+            for (int i = 0; (i < n) && issymmetric; i++) 
+            {
+                issymmetric = (A(i,j) == A(j,i));
+            }
+        }
+
+        if (issymmetric) 
+        {
+            V = A;
+
+            // Tridiagonalize.
+            tred2();
+
+            // Diagonalize.
+            tql2();
+
+        } 
+        else 
+        {
+            H = A;
+            ort.set_size(n);
+
+            // Reduce to Hessenberg form.
+            orthes();
+
+            // Reduce Hessenberg to real Schur form.
+            hqr2();
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename eigenvalue_decomposition<matrix_exp_type>::matrix_type& eigenvalue_decomposition<matrix_exp_type>::
+    get_pseudo_v (
+    ) const
+    {
+        return V;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    long eigenvalue_decomposition<matrix_exp_type>::
+    dim (
+    ) const
+    {
+        return V.nr();
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename eigenvalue_decomposition<matrix_exp_type>::complex_column_vector_type eigenvalue_decomposition<matrix_exp_type>::
+    get_eigenvalues (
+    ) const
+    {
+        return complex_matrix(get_real_eigenvalues(), get_imag_eigenvalues());
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename eigenvalue_decomposition<matrix_exp_type>::column_vector_type& eigenvalue_decomposition<matrix_exp_type>::
+    get_real_eigenvalues (
+    ) const
+    {
+        return d;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename eigenvalue_decomposition<matrix_exp_type>::column_vector_type& eigenvalue_decomposition<matrix_exp_type>::
+    get_imag_eigenvalues (
+    ) const
+    {
+        return e;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename eigenvalue_decomposition<matrix_exp_type>::complex_matrix_type eigenvalue_decomposition<matrix_exp_type>::
+    get_d (
+    ) const 
+    {
+        return diagm(complex_matrix(get_real_eigenvalues(), get_imag_eigenvalues()));
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename eigenvalue_decomposition<matrix_exp_type>::complex_matrix_type eigenvalue_decomposition<matrix_exp_type>::
+    get_v (
+    ) const 
+    {
+        complex_matrix_type CV(n,n);
+
+        for (int i = 0; i < n; i++) 
+        {
+            if (e(i) > 0) 
+            {
+                set_colm(CV,i) = complex_matrix(colm(V,i), colm(V,i+1));
+            } 
+            else if (e(i) < 0) 
+            {
+                set_colm(CV,i) = complex_matrix(colm(V,i), colm(V,i-1));
+            }
+            else
+            {
+                set_colm(CV,i) = complex_matrix(colm(V,i), uniform_matrix<type>(n,1,0));
+            }
+        }
+
+        return CV;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename eigenvalue_decomposition<matrix_exp_type>::matrix_type eigenvalue_decomposition<matrix_exp_type>::
+    get_pseudo_d (
+    ) const 
+    {
+        matrix_type D(n,n);
+
+        for (int i = 0; i < n; i++) 
+        {
+            for (int j = 0; j < n; j++) 
+            {
+                D(i,j) = 0.0;
+            }
+            D(i,i) = d(i);
+            if (e(i) > 0) 
+            {
+                D(i,i+1) = e(i);
+            } 
+            else if (e(i) < 0) 
+            {
+                D(i,i-1) = e(i);
+            }
+        }
+
+        return D;
+    }
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+//                        Private member functions
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+// Symmetric Householder reduction to tridiagonal form.
+    template <typename matrix_exp_type>
+    void eigenvalue_decomposition<matrix_exp_type>::
+    tred2() 
+    {
+        using std::abs;
+        using std::sqrt;
+
+        //  This is derived from the Algol procedures tred2 by
+        //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+        //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+        //  Fortran subroutine in EISPACK.
+
+        for (int j = 0; j < n; j++) 
+        {
+            d(j) = V(n-1,j);
+        }
+
+        // Householder reduction to tridiagonal form.
+
+        for (int i = n-1; i > 0; i--) 
+        {
+
+            // Scale to avoid under/overflow.
+
+            type scale = 0.0;
+            type h = 0.0;
+            for (int k = 0; k < i; k++) 
+            {
+                scale = scale + abs(d(k));
+            }
+            if (scale == 0.0) 
+            {
+                e(i) = d(i-1);
+                for (int j = 0; j < i; j++) 
+                {
+                    d(j) = V(i-1,j);
+                    V(i,j) = 0.0;
+                    V(j,i) = 0.0;
+                }
+            }
+            else 
+            {
+
+                // Generate Householder vector.
+
+                for (int k = 0; k < i; k++) 
+                {
+                    d(k) /= scale;
+                    h += d(k) * d(k);
+                }
+                type f = d(i-1);
+                type g = sqrt(h);
+                if (f > 0) 
+                {
+                    g = -g;
+                }
+                e(i) = scale * g;
+                h = h - f * g;
+                d(i-1) = f - g;
+                for (int j = 0; j < i; j++) 
+                {
+                    e(j) = 0.0;
+                }
+
+                // Apply similarity transformation to remaining columns.
+
+                for (int j = 0; j < i; j++) 
+                {
+                    f = d(j);
+                    V(j,i) = f;
+                    g = e(j) + V(j,j) * f;
+                    for (int k = j+1; k <= i-1; k++) 
+                    {
+                        g += V(k,j) * d(k);
+                        e(k) += V(k,j) * f;
+                    }
+                    e(j) = g;
+                }
+                f = 0.0;
+                for (int j = 0; j < i; j++) 
+                {
+                    e(j) /= h;
+                    f += e(j) * d(j);
+                }
+                type hh = f / (h + h);
+                for (int j = 0; j < i; j++) 
+                {
+                    e(j) -= hh * d(j);
+                }
+                for (int j = 0; j < i; j++) 
+                {
+                    f = d(j);
+                    g = e(j);
+                    for (int k = j; k <= i-1; k++) 
+                    {
+                        V(k,j) -= (f * e(k) + g * d(k));
+                    }
+                    d(j) = V(i-1,j);
+                    V(i,j) = 0.0;
+                }
+            }
+            d(i) = h;
+        }
+
+        // Accumulate transformations.
+
+        for (int i = 0; i < n-1; i++) 
+        {
+            V(n-1,i) = V(i,i);
+            V(i,i) = 1.0;
+            type h = d(i+1);
+            if (h != 0.0) 
+            {
+                for (int k = 0; k <= i; k++) 
+                {
+                    d(k) = V(k,i+1) / h;
+                }
+                for (int j = 0; j <= i; j++) 
+                {
+                    type g = 0.0;
+                    for (int k = 0; k <= i; k++) 
+                    {
+                        g += V(k,i+1) * V(k,j);
+                    }
+                    for (int k = 0; k <= i; k++) 
+                    {
+                        V(k,j) -= g * d(k);
+                    }
+                }
+            }
+            for (int k = 0; k <= i; k++) 
+            {
+                V(k,i+1) = 0.0;
+            }
+        }
+        for (int j = 0; j < n; j++) 
+        {
+            d(j) = V(n-1,j);
+            V(n-1,j) = 0.0;
+        }
+        V(n-1,n-1) = 1.0;
+        e(0) = 0.0;
+    } 
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    void eigenvalue_decomposition<matrix_exp_type>::
+    tql2 () 
+    {
+        using std::pow;
+        using std::min;
+        using std::max;
+        using std::abs;
+
+        //  This is derived from the Algol procedures tql2, by
+        //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+        //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+        //  Fortran subroutine in EISPACK.
+
+        for (int i = 1; i < n; i++) 
+        {
+            e(i-1) = e(i);
+        }
+        e(n-1) = 0.0;
+
+        type f = 0.0;
+        type tst1 = 0.0;
+        type eps = pow(2.0,-52.0);
+        for (int l = 0; l < n; l++) 
+        {
+
+            // Find small subdiagonal element
+
+            tst1 = max(tst1,abs(d(l)) + abs(e(l)));
+            int m = l;
+
+            // Original while-loop from Java code
+            while (m < n) 
+            {
+                if (abs(e(m)) <= eps*tst1) 
+                {
+                    break;
+                }
+                m++;
+            }
+
+
+            // If m == l, d(l) is an eigenvalue,
+            // otherwise, iterate.
+
+            if (m > l) 
+            {
+                int iter = 0;
+                do 
+                {
+                    iter = iter + 1;  // (Could check iteration count here.)
+
+                    // Compute implicit shift
+
+                    type g = d(l);
+                    type p = (d(l+1) - g) / (2.0 * e(l));
+                    type r = hypot(p,1.0);
+                    if (p < 0) 
+                    {
+                        r = -r;
+                    }
+                    d(l) = e(l) / (p + r);
+                    d(l+1) = e(l) * (p + r);
+                    type dl1 = d(l+1);
+                    type h = g - d(l);
+                    for (int i = l+2; i < n; i++) 
+                    {
+                        d(i) -= h;
+                    }
+                    f = f + h;
+
+                    // Implicit QL transformation.
+
+                    p = d(m);
+                    type c = 1.0;
+                    type c2 = c;
+                    type c3 = c;
+                    type el1 = e(l+1);
+                    type s = 0.0;
+                    type s2 = 0.0;
+                    for (int i = m-1; i >= l; i--) 
+                    {
+                        c3 = c2;
+                        c2 = c;
+                        s2 = s;
+                        g = c * e(i);
+                        h = c * p;
+                        r = hypot(p,e(i));
+                        e(i+1) = s * r;
+                        s = e(i) / r;
+                        c = p / r;
+                        p = c * d(i) - s * g;
+                        d(i+1) = h + s * (c * g + s * d(i));
+
+                        // Accumulate transformation.
+
+                        for (int k = 0; k < n; k++) 
+                        {
+                            h = V(k,i+1);
+                            V(k,i+1) = s * V(k,i) + c * h;
+                            V(k,i) = c * V(k,i) - s * h;
+                        }
+                    }
+                    p = -s * s2 * c3 * el1 * e(l) / dl1;
+                    e(l) = s * p;
+                    d(l) = c * p;
+
+                    // Check for convergence.
+
+                } while (abs(e(l)) > eps*tst1);
+            }
+            d(l) = d(l) + f;
+            e(l) = 0.0;
+        }
+
+        /*
+            The code to sort the eigenvalues and eigenvectors 
+            has been removed from here since, in the non-symmetric case,
+            we can't sort the eigenvalues in a meaningful way.  If we left this
+            code in here then the user might supply what they thought was a symmetric
+            matrix but was actually slightly non-symmetric due to rounding error
+            and then they would end up in the non-symmetric eigenvalue solver
+            where the eigenvalues don't end up getting sorted.  So to avoid
+            any possible user confusion I'm just removing this.
+        */
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    void eigenvalue_decomposition<matrix_exp_type>::
+    orthes () 
+    {
+        using std::abs;
+        using std::sqrt;
+
+        //  This is derived from the Algol procedures orthes and ortran,
+        //  by Martin and Wilkinson, Handbook for Auto. Comp.,
+        //  Vol.ii-Linear Algebra, and the corresponding
+        //  Fortran subroutines in EISPACK.
+
+        int low = 0;
+        int high = n-1;
+
+        for (int m = low+1; m <= high-1; m++) 
+        {
+
+            // Scale column.
+
+            type scale = 0.0;
+            for (int i = m; i <= high; i++) 
+            {
+                scale = scale + abs(H(i,m-1));
+            }
+            if (scale != 0.0) 
+            {
+
+                // Compute Householder transformation.
+
+                type h = 0.0;
+                for (int i = high; i >= m; i--) 
+                {
+                    ort(i) = H(i,m-1)/scale;
+                    h += ort(i) * ort(i);
+                }
+                type g = sqrt(h);
+                if (ort(m) > 0) 
+                {
+                    g = -g;
+                }
+                h = h - ort(m) * g;
+                ort(m) = ort(m) - g;
+
+                // Apply Householder similarity transformation
+                // H = (I-u*u'/h)*H*(I-u*u')/h)
+
+                for (int j = m; j < n; j++) 
+                {
+                    type f = 0.0;
+                    for (int i = high; i >= m; i--) 
+                    {
+                        f += ort(i)*H(i,j);
+                    }
+                    f = f/h;
+                    for (int i = m; i <= high; i++) 
+                    {
+                        H(i,j) -= f*ort(i);
+                    }
+                }
+
+                for (int i = 0; i <= high; i++) 
+                {
+                    type f = 0.0;
+                    for (int j = high; j >= m; j--) 
+                    {
+                        f += ort(j)*H(i,j);
+                    }
+                    f = f/h;
+                    for (int j = m; j <= high; j++) 
+                    {
+                        H(i,j) -= f*ort(j);
+                    }
+                }
+                ort(m) = scale*ort(m);
+                H(m,m-1) = scale*g;
+            }
+        }
+
+        // Accumulate transformations (Algol's ortran).
+
+        for (int i = 0; i < n; i++) 
+        {
+            for (int j = 0; j < n; j++) 
+            {
+                V(i,j) = (i == j ? 1.0 : 0.0);
+            }
+        }
+
+        for (int m = high-1; m >= low+1; m--) 
+        {
+            if (H(m,m-1) != 0.0) 
+            {
+                for (int i = m+1; i <= high; i++) 
+                {
+                    ort(i) = H(i,m-1);
+                }
+                for (int j = m; j <= high; j++) 
+                {
+                    type g = 0.0;
+                    for (int i = m; i <= high; i++) 
+                    {
+                        g += ort(i) * V(i,j);
+                    }
+                    // Double division avoids possible underflow
+                    g = (g / ort(m)) / H(m,m-1);
+                    for (int i = m; i <= high; i++) 
+                    {
+                        V(i,j) += g * ort(i);
+                    }
+                }
+            }
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    void eigenvalue_decomposition<matrix_exp_type>::
+    cdiv(type xr, type xi, type yr, type yi)  
+    {
+        using std::abs;
+        type r,d;
+        if (abs(yr) > abs(yi)) 
+        {
+            r = yi/yr;
+            d = yr + r*yi;
+            cdivr = (xr + r*xi)/d;
+            cdivi = (xi - r*xr)/d;
+        } 
+        else 
+        {
+            r = yr/yi;
+            d = yi + r*yr;
+            cdivr = (r*xr + xi)/d;
+            cdivi = (r*xi - xr)/d;
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    void eigenvalue_decomposition<matrix_exp_type>::
+    hqr2 ()
+    {
+        using std::pow;
+        using std::min;
+        using std::max;
+        using std::abs;
+        using std::sqrt;
+
+        //  This is derived from the Algol procedure hqr2,
+        //  by Martin and Wilkinson, Handbook for Auto. Comp.,
+        //  Vol.ii-Linear Algebra, and the corresponding
+        //  Fortran subroutine in EISPACK.
+
+        // Initialize
+
+        int nn = this->n;
+        int n = nn-1;
+        int low = 0;
+        int high = nn-1;
+        type eps = pow(2.0,-52.0);
+        type exshift = 0.0;
+        type p=0,q=0,r=0,s=0,z=0,t,w,x,y;
+
+        // Store roots isolated by balanc and compute matrix norm
+
+        type norm = 0.0;
+        for (int i = 0; i < nn; i++) 
+        {
+            if ((i < low) || (i > high)) 
+            {
+                d(i) = H(i,i);
+                e(i) = 0.0;
+            }
+            for (int j = max(i-1,0); j < nn; j++) 
+            {
+                norm = norm + abs(H(i,j));
+            }
+        }
+
+        // Outer loop over eigenvalue index
+
+        int iter = 0;
+        while (n >= low) 
+        {
+
+            // Look for single small sub-diagonal element
+
+            int l = n;
+            while (l > low) 
+            {
+                s = abs(H(l-1,l-1)) + abs(H(l,l));
+                if (s == 0.0) 
+                {
+                    s = norm;
+                }
+                if (abs(H(l,l-1)) < eps * s) 
+                {
+                    break;
+                }
+                l--;
+            }
+
+            // Check for convergence
+            // One root found
+
+            if (l == n) 
+            {
+                H(n,n) = H(n,n) + exshift;
+                d(n) = H(n,n);
+                e(n) = 0.0;
+                n--;
+                iter = 0;
+
+                // Two roots found
+
+            } 
+            else if (l == n-1) 
+            {
+                w = H(n,n-1) * H(n-1,n);
+                p = (H(n-1,n-1) - H(n,n)) / 2.0;
+                q = p * p + w;
+                z = sqrt(abs(q));
+                H(n,n) = H(n,n) + exshift;
+                H(n-1,n-1) = H(n-1,n-1) + exshift;
+                x = H(n,n);
+
+                // type pair
+
+                if (q >= 0) 
+                {
+                    if (p >= 0) 
+                    {
+                        z = p + z;
+                    } 
+                    else 
+                    {
+                        z = p - z;
+                    }
+                    d(n-1) = x + z;
+                    d(n) = d(n-1);
+                    if (z != 0.0) 
+                    {
+                        d(n) = x - w / z;
+                    }
+                    e(n-1) = 0.0;
+                    e(n) = 0.0;
+                    x = H(n,n-1);
+                    s = abs(x) + abs(z);
+                    p = x / s;
+                    q = z / s;
+                    r = sqrt(p * p+q * q);
+                    p = p / r;
+                    q = q / r;
+
+                    // Row modification
+
+                    for (int j = n-1; j < nn; j++) 
+                    {
+                        z = H(n-1,j);
+                        H(n-1,j) = q * z + p * H(n,j);
+                        H(n,j) = q * H(n,j) - p * z;
+                    }
+
+                    // Column modification
+
+                    for (int i = 0; i <= n; i++) 
+                    {
+                        z = H(i,n-1);
+                        H(i,n-1) = q * z + p * H(i,n);
+                        H(i,n) = q * H(i,n) - p * z;
+                    }
+
+                    // Accumulate transformations
+
+                    for (int i = low; i <= high; i++) 
+                    {
+                        z = V(i,n-1);
+                        V(i,n-1) = q * z + p * V(i,n);
+                        V(i,n) = q * V(i,n) - p * z;
+                    }
+
+                    // Complex pair
+
+                } 
+                else 
+                {
+                    d(n-1) = x + p;
+                    d(n) = x + p;
+                    e(n-1) = z;
+                    e(n) = -z;
+                }
+                n = n - 2;
+                iter = 0;
+
+                // No convergence yet
+
+            } 
+            else 
+            {
+
+                // Form shift
+
+                x = H(n,n);
+                y = 0.0;
+                w = 0.0;
+                if (l < n) 
+                {
+                    y = H(n-1,n-1);
+                    w = H(n,n-1) * H(n-1,n);
+                }
+
+                // Wilkinson's original ad hoc shift
+
+                if (iter == 10) 
+                {
+                    exshift += x;
+                    for (int i = low; i <= n; i++) 
+                    {
+                        H(i,i) -= x;
+                    }
+                    s = abs(H(n,n-1)) + abs(H(n-1,n-2));
+                    x = y = 0.75 * s;
+                    w = -0.4375 * s * s;
+                }
+
+                // MATLAB's new ad hoc shift
+
+                if (iter == 30) 
+                {
+                    s = (y - x) / 2.0;
+                    s = s * s + w;
+                    if (s > 0) 
+                    {
+                        s = sqrt(s);
+                        if (y < x) 
+                        {
+                            s = -s;
+                        }
+                        s = x - w / ((y - x) / 2.0 + s);
+                        for (int i = low; i <= n; i++) 
+                        {
+                            H(i,i) -= s;
+                        }
+                        exshift += s;
+                        x = y = w = 0.964;
+                    }
+                }
+
+                iter = iter + 1;   // (Could check iteration count here.)
+
+                // Look for two consecutive small sub-diagonal elements
+
+                int m = n-2;
+                while (m >= l) 
+                {
+                    z = H(m,m);
+                    r = x - z;
+                    s = y - z;
+                    p = (r * s - w) / H(m+1,m) + H(m,m+1);
+                    q = H(m+1,m+1) - z - r - s;
+                    r = H(m+2,m+1);
+                    s = abs(p) + abs(q) + abs(r);
+                    p = p / s;
+                    q = q / s;
+                    r = r / s;
+                    if (m == l) 
+                    {
+                        break;
+                    }
+                    if (abs(H(m,m-1)) * (abs(q) + abs(r)) <
+                        eps * (abs(p) * (abs(H(m-1,m-1)) + abs(z) +
+                                         abs(H(m+1,m+1))))) 
+                    {
+                        break;
+                    }
+                    m--;
+                }
+
+                for (int i = m+2; i <= n; i++) 
+                {
+                    H(i,i-2) = 0.0;
+                    if (i > m+2) 
+                    {
+                        H(i,i-3) = 0.0;
+                    }
+                }
+
+                // Double QR step involving rows l:n and columns m:n
+
+                for (int k = m; k <= n-1; k++) 
+                {
+                    int notlast = (k != n-1);
+                    if (k != m) 
+                    {
+                        p = H(k,k-1);
+                        q = H(k+1,k-1);
+                        r = (notlast ? H(k+2,k-1) : 0.0);
+                        x = abs(p) + abs(q) + abs(r);
+                        if (x != 0.0) 
+                        {
+                            p = p / x;
+                            q = q / x;
+                            r = r / x;
+                        }
+                    }
+                    if (x == 0.0) 
+                    {
+                        break;
+                    }
+                    s = sqrt(p * p + q * q + r * r);
+                    if (p < 0) 
+                    {
+                        s = -s;
+                    }
+                    if (s != 0) 
+                    {
+                        if (k != m) 
+                        {
+                            H(k,k-1) = -s * x;
+                        } 
+                        else if (l != m) 
+                        {
+                            H(k,k-1) = -H(k,k-1);
+                        }
+                        p = p + s;
+                        x = p / s;
+                        y = q / s;
+                        z = r / s;
+                        q = q / p;
+                        r = r / p;
+
+                        // Row modification
+
+                        for (int j = k; j < nn; j++) 
+                        {
+                            p = H(k,j) + q * H(k+1,j);
+                            if (notlast) 
+                            {
+                                p = p + r * H(k+2,j);
+                                H(k+2,j) = H(k+2,j) - p * z;
+                            }
+                            H(k,j) = H(k,j) - p * x;
+                            H(k+1,j) = H(k+1,j) - p * y;
+                        }
+
+                        // Column modification
+
+                        for (int i = 0; i <= min(n,k+3); i++) 
+                        {
+                            p = x * H(i,k) + y * H(i,k+1);
+                            if (notlast) 
+                            {
+                                p = p + z * H(i,k+2);
+                                H(i,k+2) = H(i,k+2) - p * r;
+                            }
+                            H(i,k) = H(i,k) - p;
+                            H(i,k+1) = H(i,k+1) - p * q;
+                        }
+
+                        // Accumulate transformations
+
+                        for (int i = low; i <= high; i++) 
+                        {
+                            p = x * V(i,k) + y * V(i,k+1);
+                            if (notlast) 
+                            {
+                                p = p + z * V(i,k+2);
+                                V(i,k+2) = V(i,k+2) - p * r;
+                            }
+                            V(i,k) = V(i,k) - p;
+                            V(i,k+1) = V(i,k+1) - p * q;
+                        }
+                    }  // (s != 0)
+                }  // k loop
+            }  // check convergence
+        }  // while (n >= low)
+
+        // Backsubstitute to find vectors of upper triangular form
+
+        if (norm == 0.0) 
+        {
+            return;
+        }
+
+        for (n = nn-1; n >= 0; n--) 
+        {
+            p = d(n);
+            q = e(n);
+
+            // Real vector
+
+            if (q == 0) 
+            {
+                int l = n;
+                H(n,n) = 1.0;
+                for (int i = n-1; i >= 0; i--) 
+                {
+                    w = H(i,i) - p;
+                    r = 0.0;
+                    for (int j = l; j <= n; j++) 
+                    {
+                        r = r + H(i,j) * H(j,n);
+                    }
+                    if (e(i) < 0.0) 
+                    {
+                        z = w;
+                        s = r;
+                    } 
+                    else 
+                    {
+                        l = i;
+                        if (e(i) == 0.0) 
+                        {
+                            if (w != 0.0) 
+                            {
+                                H(i,n) = -r / w;
+                            } 
+                            else 
+                            {
+                                H(i,n) = -r / (eps * norm);
+                            }
+
+                            // Solve real equations
+
+                        } 
+                        else 
+                        {
+                            x = H(i,i+1);
+                            y = H(i+1,i);
+                            q = (d(i) - p) * (d(i) - p) + e(i) * e(i);
+                            t = (x * s - z * r) / q;
+                            H(i,n) = t;
+                            if (abs(x) > abs(z)) 
+                            {
+                                H(i+1,n) = (-r - w * t) / x;
+                            } 
+                            else 
+                            {
+                                H(i+1,n) = (-s - y * t) / z;
+                            }
+                        }
+
+                        // Overflow control
+
+                        t = abs(H(i,n));
+                        if ((eps * t) * t > 1) 
+                        {
+                            for (int j = i; j <= n; j++) 
+                            {
+                                H(j,n) = H(j,n) / t;
+                            }
+                        }
+                    }
+                }
+
+                // Complex vector
+
+            } 
+            else if (q < 0) 
+            {
+                int l = n-1;
+
+                // Last vector component imaginary so matrix is triangular
+
+                if (abs(H(n,n-1)) > abs(H(n-1,n))) 
+                {
+                    H(n-1,n-1) = q / H(n,n-1);
+                    H(n-1,n) = -(H(n,n) - p) / H(n,n-1);
+                } 
+                else 
+                {
+                    cdiv(0.0,-H(n-1,n),H(n-1,n-1)-p,q);
+                    H(n-1,n-1) = cdivr;
+                    H(n-1,n) = cdivi;
+                }
+                H(n,n-1) = 0.0;
+                H(n,n) = 1.0;
+                for (int i = n-2; i >= 0; i--) 
+                {
+                    type ra,sa,vr,vi;
+                    ra = 0.0;
+                    sa = 0.0;
+                    for (int j = l; j <= n; j++) 
+                    {
+                        ra = ra + H(i,j) * H(j,n-1);
+                        sa = sa + H(i,j) * H(j,n);
+                    }
+                    w = H(i,i) - p;
+
+                    if (e(i) < 0.0) 
+                    {
+                        z = w;
+                        r = ra;
+                        s = sa;
+                    } 
+                    else 
+                    {
+                        l = i;
+                        if (e(i) == 0) 
+                        {
+                            cdiv(-ra,-sa,w,q);
+                            H(i,n-1) = cdivr;
+                            H(i,n) = cdivi;
+                        } 
+                        else 
+                        {
+
+                            // Solve complex equations
+
+                            x = H(i,i+1);
+                            y = H(i+1,i);
+                            vr = (d(i) - p) * (d(i) - p) + e(i) * e(i) - q * q;
+                            vi = (d(i) - p) * 2.0 * q;
+                            if ((vr == 0.0) && (vi == 0.0)) 
+                            {
+                                vr = eps * norm * (abs(w) + abs(q) +
+                                                   abs(x) + abs(y) + abs(z));
+                            }
+                            cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
+                            H(i,n-1) = cdivr;
+                            H(i,n) = cdivi;
+                            if (abs(x) > (abs(z) + abs(q))) 
+                            {
+                                H(i+1,n-1) = (-ra - w * H(i,n-1) + q * H(i,n)) / x;
+                                H(i+1,n) = (-sa - w * H(i,n) - q * H(i,n-1)) / x;
+                            }
+                            else 
+                            {
+                                cdiv(-r-y*H(i,n-1),-s-y*H(i,n),z,q);
+                                H(i+1,n-1) = cdivr;
+                                H(i+1,n) = cdivi;
+                            }
+                        }
+
+                        // Overflow control
+
+                        t = max(abs(H(i,n-1)),abs(H(i,n)));
+                        if ((eps * t) * t > 1) 
+                        {
+                            for (int j = i; j <= n; j++) 
+                            {
+                                H(j,n-1) = H(j,n-1) / t;
+                                H(j,n) = H(j,n) / t;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        // Vectors of isolated roots
+
+        for (int i = 0; i < nn; i++) 
+        {
+            if (i < low || i > high) 
+            {
+                for (int j = i; j < nn; j++) 
+                {
+                    V(i,j) = H(i,j);
+                }
+            }
+        }
+
+        // Back transformation to get eigenvectors of original matrix
+
+        for (int j = nn-1; j >= low; j--) 
+        {
+            for (int i = low; i <= high; i++) 
+            {
+                z = 0.0;
+                for (int k = low; k <= min(j,high); k++) 
+                {
+                    z = z + V(i,k) * H(k,j);
+                }
+                V(i,j) = z;
+            }
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+
+} 
+
+#endif // DLIB_MATRIX_EIGENVALUE_DECOMPOSITION_H 
+
+
+
+

dlib/matrix/matrix_la.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_MATRIx_LA_FUNCTS_
+#define DLIB_MATRIx_LA_FUNCTS_ 
+
+#include "matrix_la_abstract.h"
+
+#include "matrix_lu.h"
+#include "matrix_qr.h"
+#include "matrix_cholesky.h"
+#include "matrix_eigenvalue.h"
+
+namespace dlib
+{
+
+}
+
+#endif // DLIB_MATRIx_LA_FUNCTS_
+
+

dlib/matrix/matrix_la_abstract.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_
+#ifdef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_ 
+
+#include "matrix_abstract.h"
+#include <complex>
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename matrix_exp_type
+        >
+    class lu_decomposition
+    {
+        /*!
+            REQUIREMENTS ON matrix_exp_type
+                must be some kind of matrix expression as defined in the 
+                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
+
+            WHAT THIS OBJECT REPRESENTS
+                This object represents something that can compute an LU 
+                decomposition of a real valued matrix.  That is, for any 
+                matrix A it computes matrices L, U, and a pivot vector P such 
+                that rowm(A,P) == L*U.
+
+                The LU decomposition with pivoting always exists, even if the matrix is
+                singular, so the constructor will never fail.  The primary use of the
+                LU decomposition is in the solution of square systems of simultaneous
+                linear equations.  This will fail if is_singular() returns true (or
+                if A is very nearly singular).
+        !*/
+
+    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;
+        typedef matrix<long,NR,1,mem_manager_type,layout_type> pivot_column_vector_type;
+
+        template <typename EXP>
+        lu_decomposition (
+            const matrix_exp<EXP> &A
+        );
+        /*!
+            requires
+                - EXP::type == lu_decomposition::type 
+                - A.size() > 0
+            ensures
+                - #nr() == A.nr()
+                - #nc() == A.nc()
+                - #is_square() == (A.nr() == A.nc())
+                - computes the LU factorization of the given A matrix.
+        !*/
+
+        bool is_square (
+        ) const;
+        /*!
+            ensures
+                - if (the input A matrix was a square matrix) then
+                    - returns true
+                - else
+                    - returns false
+        !*/
+
+        bool is_singular (
+        ) const;
+        /*!
+            requires
+                - is_square() == true
+            ensures
+                - if (the input A matrix is singular) then
+                    - returns true
+                - else
+                    - returns false
+        !*/
+
+        long nr(
+        ) const;
+        /*!
+            ensures
+                - returns the number of rows in the input matrix
+        !*/
+
+        long nc(
+        ) const;
+        /*!
+            ensures
+                - returns the number of columns in the input matrix
+        !*/
+
+        const matrix_type get_l (
+        ) const; 
+        /*!
+            ensures
+                - returns the lower triangular L factor of the LU factorization.  
+                - L.nr() == nr()
+                - L.nc() == min(nr(),nc())
+        !*/
+
+        const matrix_type get_u (
+        ) const;
+        /*!
+            ensures
+                - returns the upper triangular U factor of the LU factorization.  
+                - U.nr() == min(nr(),nc())
+                - U.nc() == nc()
+        !*/
+
+        const pivot_column_vector_type& get_pivot (
+        ) const;
+        /*!
+            ensures
+                - returns the pivot permutation vector.  That is,
+                  if A is the input matrix then this function 
+                  returns a vector P such that:
+                    - rowm(A,P) == get_l()*get_u() 
+                    - P.nr() == A.nr()
+        !*/
+
+        type det (
+        ) const;
+        /*!
+            requires
+                - is_square() == true
+            ensures
+                - computes and returns the determinant of the input 
+                  matrix using LU factors.
+        !*/
+
+        template <typename EXP>
+        const matrix_type solve (
+            const matrix_exp<EXP> &B
+        ) const;
+        /*!
+            requires
+                - EXP::type == lu_decomposition::type
+                - is_square() == true
+                - B.nr() == nr()
+            ensures
+                - Let A denote the input matrix to this class's constructor.  
+                  Then this function solves A*X == B for X and returns X.  
+                - Note that if A is singular (or very close to singular) then
+                  the X returned by this function won't fit A*X == B very well (if at all).
+        !*/
+
+    }; 
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename matrix_exp_type
+        >
+    class cholesky_decomposition
+    {
+        /*! 
+            REQUIREMENTS ON matrix_exp_type
+                must be some kind of matrix expression as defined in the 
+                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
+
+            WHAT THIS OBJECT REPRESENTS
+                This object represents something that can compute a cholesky 
+                decomposition of a real valued matrix.  That is, for any 
+                symmetric, positive definite matrix A, it computes a lower 
+                triangular matrix L such that A == L*trans(L).
+                
+                If the matrix is not symmetric or positive definite, the function
+                computes only a partial decomposition.  This can be tested with
+                the is_spd() flag.
+        !*/
+
+    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 typename matrix_exp_type::matrix_type matrix_type;
+        typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;
+
+        template <typename EXP>
+        cholesky_decomposition(
+            const matrix_exp<EXP>& A
+        );
+        /*!
+            requires
+                - EXP::type == cholesky_decomposition::type 
+                - A.size() > 0
+                - A.nr() == A.nc() 
+                  (i.e. A must be a square matrix)
+            ensures
+                - if (A is symmetric positive-definite) then
+                    - #is_spd() == true 
+                    - Constructs a lower triangular matrix L, such that L*trans(L) == A.
+                      and #get_l() == L
+                - else
+                    - #is_spd() == false
+        !*/
+
+        bool is_spd(
+        ) const;
+        /*!
+            ensures
+                - if (the input matrix was symmetric positive-definite) then
+                    - returns true
+                - else
+                    - returns false
+        !*/
+
+        const matrix_type& get_l(
+        ) const;
+        /*!
+            ensures
+                - returns the lower triangular factor, L, such that L*trans(L) == A
+                  (where A is the input matrix to this class's constructor)
+                - Note that if A is not symmetric positive definite or positive semi-definite
+                  then the equation L*trans(L) == A won't hold.  
+        !*/
+
+        template <typename EXP>
+        const matrix solve (
+            const matrix_exp<EXP>& B
+        ) const;
+        /*!
+            requires
+                - EXP::type == cholesky_decomposition::type
+                - B.nr() == get_l().nr()
+                  (i.e. the number of rows in B must match the number of rows in the
+                  input matrix A)
+            ensures
+                - Let A denote the input matrix to this class's constructor.  Then 
+                  this function solves A*X = B for X and returns X.  
+                - Note that if is_spd() == false or A was really close to being
+                  non-SPD then the solver will fail to find an accurate solution.
+        !*/
+
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename matrix_exp_type
+        >
+    class qr_decomposition 
+    {
+        /*! 
+            REQUIREMENTS ON matrix_exp_type
+                must be some kind of matrix expression as defined in the 
+                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
+
+            WHAT THIS OBJECT REPRESENTS
+                This object represents something that can compute a classical
+                QR decomposition of an m-by-n real valued matrix A with m >= n.  
+
+                The QR decomposition is an m-by-n orthogonal matrix Q and an 
+                n-by-n upper triangular matrix R so that A == Q*R. The QR decomposition 
+                always exists, even if the matrix does not have full rank, so the 
+                constructor will never fail.  The primary use of the QR decomposition 
+                is in the least squares solution of non-square systems of simultaneous 
+                linear equations.  This will fail if is_full_rank() returns false or
+                A is very nearly not full rank.
+
+                The Q and R factors can be retrieved via the get_q() and get_r()
+                methods. Furthermore, a solve() method is provided to find the
+                least squares solution of Ax=b using the QR factors.  
+        !*/
+
+    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,0,1,mem_manager_type,layout_type> column_vector_type;
+
+        template <typename EXP>
+        qr_decomposition(
+            const matrix_exp<EXP>& A
+        );
+        /*!
+            requires
+                - EXP::type == qr_decomposition::type
+                - A.nr() >= A.nc()
+                - A.size() > 0
+            ensures
+                - #nr() == A.nr()
+                - #nc() == A.nc()
+                - computes the QR decomposition of the given A matrix.
+        !*/
+
+        bool is_full_rank(
+        ) const;
+        /*!
+            ensures
+                - if (the input A matrix had full rank) then
+                    - returns true
+                - else
+                    - returns false
+        !*/
+
+        long nr(
+        ) const;
+        /*!
+            ensures
+                - returns the number of rows in the input matrix
+        !*/
+
+        long nc(
+        ) const;
+        /*!
+            ensures
+                - returns the number of columns in the input matrix
+        !*/
+
+        const matrix_type get_householder (
+        )  const;
+        /*!
+            ensures
+                - returns a matrix H such that:
+                    - H is the lower trapezoidal matrix whose columns define the 
+                      Householder reflection vectors from QR factorization 
+                    - H.nr() == nr()
+                    - H.nc() == nc()
+        !*/
+
+        const matrix_type get_r (
+        ) const;
+        /*!
+            ensures
+                - returns a matrix R such that: 
+                    - R is the upper triangular factor, R, of the QR factorization
+                    - get_q()*R == input matrix A
+                    - R.nr() == nc()
+                    - R.nc() == nc()
+        !*/
+
+        const matrix_type get_q (
+        ) const;
+        /*!
+            ensures
+                - returns a matrix Q such that:
+                    - Q is the economy-sized orthogonal factor Q from the QR 
+                      factorization.  
+                    - trans(Q)*Q == identity matrix
+                    - Q*get_r() == input matrix A 
+                    - Q.nr() == nr()
+                    - Q.nc() == nc()
+        !*/
+
+        template <typename EXP>
+        const matrix_type solve (
+            const matrix_exp<EXP>& B
+        ) const;
+        /*!
+            requires
+                - EXP::type == qr_decomposition::type
+                - B.nr() == nr()
+            ensures
+                - Let A denote the input matrix to this class's constructor.  
+                  Then this function finds the least squares solution to the equation A*X = B 
+                  and returns X.  X has the following properties: 
+                    - X is the matrix that minimizes the two norm of A*X-B.  That is, it
+                      minimizes sum(squared(A*X - B)).
+                    - X.nr() == nc()
+                    - X.nc() == B.nc()
+                - Note that this function will fail to output a good solution if is_full_rank() == false
+                  or the A matrix is close to not being full rank.
+        !*/
+
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename matrix_exp_type
+        >
+    class eigenvalue_decomposition
+    {
+        /*!
+            REQUIREMENTS ON matrix_exp_type
+                must be some kind of matrix expression as defined in the 
+                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
+
+            WHAT THIS OBJECT REPRESENTS
+                This object represents something that can compute an eigenvalue 
+                decomposition of a real valued matrix.   So it gives 
+                you the set of eigenvalues and eigenvectors for a matrix.   
+
+                Let A denote the input matrix to this object's constructor.  Then 
+                what this object does is it finds two matrices, D and V, such that
+                    - A*V == V*D
+                Where V is a square matrix that contains all the eigenvectors
+                of the A matrix (each column of V is an eigenvector) and
+                D is a diagonal matrix containing the eigenvalues of A.
+
+
+                It is important to note that if A is symmetric or non-symmetric you
+                get somewhat different results.  If A is a symmetric matrix (i.e. A == trans(A))
+                then:
+                    - All the eigenvalues and eigenvectors of A are real numbers. 
+                        - Because of this there isn't really any point in using the
+                          part of this class's interface that returns complex matrices.
+                          All you need are the get_real_eigenvalues() and
+                          get_pseudo_v() functions.  
+                    - V*trans(V) should be equal to the identity matrix.  That is, all the
+                      eigenvectors in V should be linearly independent. 
+                        - So A == V*D*trans(V)
+
+                On the other hand, if A is not symmetric then:
+                    - Some of the eigenvalues and eigenvectors might be complex numbers.  
+                        - An eigenvalue is complex if and only if its corresponding eigenvector 
+                          is complex.  So you can check for this case by just checking 
+                          get_imag_eigenvalues() to see if any values are non-zero.  You don't 
+                          have to check the V matrix as well.
+                    - V*trans(V) won't be equal to the identity matrix but it is usually
+                      invertible.  So A == V*D*inv(V) is usually a valid statement but
+                      A == V*D*trans(V) won't be.
+        !*/
+
+    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 typename matrix_exp_type::matrix_type matrix_type;
+        typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;
+
+        typedef matrix<std::complex<type>,0,0,mem_manager_type,layout_type> complex_matrix_type;
+        typedef matrix<std::complex<type>,NR,1,mem_manager_type,layout_type> complex_column_vector_type;
+
+
+        template <typename EXP>
+        eigenvalue_decomposition(
+            const matrix_exp<EXP>& A
+        ); 
+        /*!
+            requires
+                - A.nr() == A.nc() 
+                - A.size() > 0
+                - EXP::type == eigenvalue_decomposition::type 
+            ensures
+                - #dim() == A.nr()
+                - computes the eigenvalue decomposition of A.  
+                - #get_eigenvalues() == the eigenvalues of A
+                - #get_v() == all the eigenvectors of A
+        !*/
+
+        long dim (
+        ) const;
+        /*!
+            ensures
+                - dim() == the number of rows/columns in the input matrix A 
+        !*/
+
+        const complex_column_vector_type get_eigenvalues (
+        ) const;
+        /*!
+            ensures
+                - returns diag(get_d()).  That is, returns a 
+                  vector that contains the eigenvalues of the input 
+                  matrix.
+                - the returned vector has dim() rows
+                - the eigenvalues are not sorted in any particular way
+        !*/
+
+        const column_vector_type& get_real_eigenvalues (
+        ) const;
+        /*! 
+            ensures
+                - returns the real parts of the eigenvalues.  That is,
+                  returns real(get_eigenvalues()) 
+                - the returned vector has dim() rows
+                - the eigenvalues are not sorted in any particular way
+        !*/
+
+        const column_vector_type& get_imag_eigenvalues (
+        ) const;
+        /*! 
+            ensures
+                - returns the imaginary parts of the eigenvalues.  That is,
+                  returns imag(get_eigenvalues()) 
+                - the returned vector has dim() rows
+                - the eigenvalues are not sorted in any particular way
+        !*/
+
+        const complex_matrix_type get_v (
+        ) const;
+        /*!
+            ensures
+                - returns the eigenvector matrix V that is 
+                  dim() rows by dim() columns
+                - Each column in V is one of the eigenvectors of the input 
+                  matrix
+        !*/
+
+        const complex_matrix_type get_d (
+        ) const; 
+        /*!
+            ensures
+                - returns a matrix D such that:
+                    - D.nr() == dim()
+                    - D.nc() == dim()
+                    - diag(D) == get_eigenvalues()
+                      (i.e. the diagonal of D contains all the eigenvalues in the input matrix)
+                    - all off diagonal elements of D are set to 0
+        !*/
+
+        const matrix_type& get_pseudo_v (
+        ) const;
+        /*!
+            ensures
+                - returns a matrix that is dim() rows by dim() columns
+                - Let A denote the input matrix given to this object's constructor.
+                - if (A has any imaginary eigenvalues) then
+                    - returns the pseudo-eigenvector matrix V  
+                    - The matrix V returned by this function is structured such that:
+                        - A*V == V*get_pseudo_d()
+                - else
+                    - returns the eigenvector matrix V with A's eigenvectors as
+                      the columns of V
+                    - A*V == V*diagm(get_real_eigenvalues())
+        !*/
+
+        const matrix_type get_pseudo_d (
+        ) const; 
+        /*!
+            ensures
+                - The returned matrix is dim() rows by dim() columns
+                - Computes and returns the block diagonal eigenvalue matrix.
+                  If the original matrix A is not symmetric, then the eigenvalue 
+                  matrix D is block diagonal with the real eigenvalues in 1-by-1 
+                  blocks and any complex eigenvalues,
+                  a + i*b, in 2-by-2 blocks, (a, b; -b, a).  That is, if the complex
+                  eigenvalues look like
+
+                      u + iv     .        .          .      .    .
+                        .      u - iv     .          .      .    .
+                        .        .      a + ib       .      .    .
+                        .        .        .        a - ib   .    .
+                        .        .        .          .      x    .
+                        .        .        .          .      .    y
+
+                  Then D looks like
+
+                        u        v        .          .      .    .
+                       -v        u        .          .      .    . 
+                        .        .        a          b      .    .
+                        .        .       -b          a      .    .
+                        .        .        .          .      x    .
+                        .        .        .          .      .    y
+
+                  This keeps V (The V you get from get_pseudo_v()) a real matrix in both 
+                  symmetric and non-symmetric cases, and A*V = V*D.
+                - the eigenvalues are not sorted in any particular way
+        !*/
+
+    };
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_MATRIx_LA_FUNCTS_ABSTRACT_
+

dlib/matrix/matrix_lu.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.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_LU_DECOMPOSITION_H
+#define DLIB_MATRIX_LU_DECOMPOSITION_H
+
+#include "matrix.h" 
+#include "matrix_utilities.h"
+#include "matrix_subexp.h"
+#include <algorithm>
+
+namespace dlib 
+{
+
+    template <
+        typename matrix_exp_type
+        >
+    class lu_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;
+        typedef matrix<long,NR,1,mem_manager_type,layout_type> pivot_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>
+        lu_decomposition (
+            const matrix_exp<EXP> &A
+        );
+
+        bool is_square (
+        ) const;
+
+        bool is_singular (
+        ) const;
+
+        long nr(
+        ) const;
+
+        long nc(
+        ) const;
+
+        const matrix_type get_l (
+        ) const; 
+
+        const matrix_type get_u (
+        ) const;
+
+        const pivot_column_vector_type& get_pivot (
+        ) const;
+
+        type det (
+        ) const;
+
+        template <typename EXP>
+        const matrix_type solve (
+            const matrix_exp<EXP> &B
+        ) const;
+
+    private:
+
+        /* Array for internal storage of decomposition.  */
+        matrix_type  LU;
+        long m, n, pivsign; 
+        pivot_column_vector_type piv;
+
+
+    }; 
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+//                              Public member functions
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    template <typename EXP>
+    lu_decomposition<matrix_exp_type>::
+    lu_decomposition (
+        const matrix_exp<EXP>& A
+    ) : 
+        LU(A),
+        m(A.nr()),
+        n(A.nc())
+    {
+        using namespace std;
+        using std::abs;
+
+        COMPILE_TIME_ASSERT((is_same_type<type, typename EXP::type>::value));
+
+        // make sure requires clause is not broken
+        DLIB_ASSERT(A.size() > 0,
+            "\tlu_decomposition::lu_decomposition(A)"
+            << "\n\tInvalid inputs were given to this function"
+            << "\n\tA.size(): " << A.size()
+            << "\n\tthis:     " << this
+            );
+
+        // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
+
+
+        piv = range(0,m-1);
+        pivsign = 1;
+
+        column_vector_type LUcolj(m);
+
+        // Outer loop.
+        for (long j = 0; j < n; j++) 
+        {
+
+            // Make a copy of the j-th column to localize references.
+            LUcolj = colm(LU,j);
+
+            // Apply previous transformations.
+            for (long i = 0; i < m; i++) 
+            {
+                // Most of the time is spent in the following dot product.
+                const long kmax = std::min(i,j);
+                const double s = rowm(LU,i, kmax)*colm(LUcolj,0,kmax);
+
+                LU(i,j) = LUcolj(i) -= s;
+            }
+
+            // Find pivot and exchange if necessary.
+            long p = j;
+            for (long i = j+1; i < m; i++) 
+            {
+                if (abs(LUcolj(i)) > abs(LUcolj(p))) 
+                {
+                    p = i;
+                }
+            }
+            if (p != j) 
+            {
+                long k=0;
+                for (k = 0; k < n; k++) 
+                {
+                    double t = LU(p,k); 
+                    LU(p,k) = LU(j,k); 
+                    LU(j,k) = t;
+                }
+                k = piv(p); 
+                piv(p) = piv(j); 
+                piv(j) = k;
+                pivsign = -pivsign;
+            }
+
+            // Compute multipliers.
+            if ((j < m) && (LU(j,j) != 0.0)) 
+            {
+                for (long i = j+1; i < m; i++) 
+                {
+                    LU(i,j) /= LU(j,j);
+                }
+            }
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    bool lu_decomposition<matrix_exp_type>::
+    is_square (
+    ) const
+    {
+        return m == n;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    long lu_decomposition<matrix_exp_type>::
+    nr (
+    ) const
+    {
+        return m;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    long lu_decomposition<matrix_exp_type>::
+    nc (
+    ) const
+    {
+        return n;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    bool lu_decomposition<matrix_exp_type>::
+    is_singular (
+    ) const
+    {
+        /* Is the matrix singular?
+          if upper triangular factor U (and hence A) is singular, false otherwise.
+        */
+        // make sure requires clause is not broken
+        DLIB_ASSERT(is_square() == true,
+            "\tbool lu_decomposition::is_singular()"
+            << "\n\tYou can only use this on square matrices"
+            << "\n\tthis: " << this
+            );
+
+        type eps = max(abs(diag(LU)));
+        if (eps != 0)
+            eps *= std::sqrt(std::numeric_limits<type>::epsilon())/10;
+        else
+            eps = 1;  // there is no max so just use 1
+
+        return min(abs(diag(LU))) < eps;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename lu_decomposition<matrix_exp_type>::matrix_type lu_decomposition<matrix_exp_type>::
+    get_l (
+    ) const
+    {
+        if (LU.nr() >= LU.nc())
+            return lowerm(LU,1.0);
+        else
+            return lowerm(subm(LU,0,0,m,m), 1.0);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename lu_decomposition<matrix_exp_type>::matrix_type lu_decomposition<matrix_exp_type>::
+    get_u (
+    ) const 
+    {
+        if (LU.nr() >= LU.nc())
+            return upperm(subm(LU,0,0,n,n));
+        else
+            return upperm(LU);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename lu_decomposition<matrix_exp_type>::pivot_column_vector_type& lu_decomposition<matrix_exp_type>::
+    get_pivot (
+    ) const
+    {
+        return piv;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    typename lu_decomposition<matrix_exp_type>::type lu_decomposition<matrix_exp_type>::
+    det (
+    ) const
+    {
+        // make sure requires clause is not broken
+        DLIB_ASSERT(is_square() == true,
+            "\ttype lu_decomposition::det()"
+            << "\n\tYou can only use this on square matrices"
+            << "\n\tthis: " << this
+            );
+
+        // Check if it is singular and if it is just return 0.  
+        // We ant to do this because a prod() operation can easily
+        // overcome a single diagonal element that is effectively 0 when
+        // LU is a big enough matrix.
+        if (is_singular())
+            return 0;
+
+        return prod(diag(LU))*static_cast<type>(pivsign);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    template <typename EXP>
+    const typename lu_decomposition<matrix_exp_type>::matrix_type lu_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(is_square() == true && B.nr() == nr(),
+            "\ttype lu_decomposition::solve()"
+            << "\n\tInvalid arguments to this function"
+            << "\n\tis_square():   " << (is_square()? "true":"false" )
+            << "\n\tB.nr():        " << B.nr() 
+            << "\n\tnr():          " << nr() 
+            << "\n\tthis:          " << this
+            );
+
+        const long nx = B.nc();
+        // if there are multiple columns in B
+        if (nx > 1)
+        {
+
+            // Copy right hand side with pivoting
+            matrix_type X(rowm(B, piv));
+
+            // Solve L*Y = B(piv,:)
+            for (long k = 0; k < n; k++) 
+            {
+                for (long i = k+1; i < n; i++) 
+                {
+                    for (long j = 0; j < nx; j++) 
+                    {
+                        X(i,j) -= X(k,j)*LU(i,k);
+                    }
+                }
+            }
+            // Solve U*X = Y;
+            for (long k = n-1; k >= 0; k--) 
+            {
+                for (long j = 0; j < nx; j++) 
+                {
+                    X(k,j) /= LU(k,k);
+                }
+                for (long i = 0; i < k; i++) 
+                {
+                    for (long j = 0; j < nx; j++) 
+                    {
+                        X(i,j) -= X(k,j)*LU(i,k);
+                    }
+                }
+            }
+            return X;
+        }
+        else
+        {
+            column_vector_type x(rowm(B, piv));
+
+            // Solve L*Y = B(piv)
+            for (long k = 0; k < n; k++) 
+            {
+                for (long i = k+1; i < n; i++) 
+                {
+                    x(i) -= x(k)*LU(i,k);
+                }
+            }
+
+            // Solve U*X = Y;
+            for (long k = n-1; k >= 0; k--) 
+            {
+                x(k) /= LU(k,k);
+                for (long i = 0; i < k; i++) 
+                    x(i) -= x(k)*LU(i,k);
+            }
+
+
+            return x;
+
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+} 
+
+#endif // DLIB_MATRIX_LU_DECOMPOSITION_H 
+
+

dlib/matrix/matrix_math_functions_abstract.h

@@ -312,6 +312,7 @@ namespace dlib
                 - R has the same dimensions as m
                 - for all valid r and c:
                   R(r,c) == std::norm(m(r,c))
+                  (note that std::norm(val) == val.real()*val.real() + val.imag()*val.imag())
     !*/
 
 // ----------------------------------------------------------------------------------------

dlib/matrix/matrix_qr.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.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_QR_DECOMPOSITION_H
+#define DLIB_MATRIX_QR_DECOMPOSITION_H
+
+#include "matrix.h" 
+#include "matrix_utilities.h"
+#include "matrix_subexp.h"
+
+namespace dlib 
+{
+
+    template <
+        typename matrix_exp_type
+        >
+    class qr_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,0,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>
+        qr_decomposition(
+            const matrix_exp<EXP>& A
+        );
+
+        bool is_full_rank(
+        ) const;
+
+        long nr(
+        ) const;
+
+        long nc(
+        ) const;
+
+        const matrix_type get_householder (
+        )  const;
+
+        const matrix_type get_r (
+        ) const;
+
+        const matrix_type get_q (
+        ) const;
+
+        template <typename EXP>
+        const matrix_type solve (
+            const matrix_exp<EXP>& B
+        ) const;
+
+    private:
+
+        template <typename EXP>
+        const matrix_type solve_mat (
+            const matrix_exp<EXP>& B
+        ) const;
+
+        template <typename EXP>
+        const matrix_type solve_vect (
+            const matrix_exp<EXP>& B
+        ) const;
+
+
+        /** Array for internal storage of decomposition.
+        @serial internal array storage.
+        */
+        matrix_type QR_;
+
+        /** Row and column dimensions.
+        @serial column dimension.
+        @serial row dimension.
+        */
+        long m, n;
+
+        /** Array for internal storage of diagonal of R.
+        @serial diagonal of R.
+        */
+        column_vector_type Rdiag;
+
+
+    };
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+//                                      Member functions
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    template <typename EXP>
+    qr_decomposition<matrix_exp_type>::
+    qr_decomposition(
+        const matrix_exp<EXP>& A
+    )		
+    {
+        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,
+            "\tqr_decomposition::qr_decomposition(A)"
+            << "\n\tInvalid inputs were given to this function"
+            << "\n\tA.nr():   " << A.nr()
+            << "\n\tA.nc():   " << A.nc()
+            << "\n\tA.size(): " << A.size()
+            << "\n\tthis:     " << this
+            );
+
+
+        QR_ = A;
+        m = A.nr();
+        n = A.nc();
+        Rdiag.set_size(n);
+        long i=0, j=0, k=0;
+
+        // Main loop.
+        for (k = 0; k < n; k++) 
+        {
+            // Compute 2-norm of k-th column without under/overflow.
+            type nrm = 0;
+            for (i = k; i < m; i++) 
+            {
+                nrm = hypot(nrm,QR_(i,k));
+            }
+
+            if (nrm != 0.0) 
+            {
+                // Form k-th Householder vector.
+                if (QR_(k,k) < 0) 
+                {
+                    nrm = -nrm;
+                }
+                for (i = k; i < m; i++) 
+                {
+                    QR_(i,k) /= nrm;
+                }
+                QR_(k,k) += 1.0;
+
+                // Apply transformation to remaining columns.
+                for (j = k+1; j < n; j++) 
+                {
+                    type s = 0.0; 
+                    for (i = k; i < m; i++) 
+                    {
+                        s += QR_(i,k)*QR_(i,j);
+                    }
+                    s = -s/QR_(k,k);
+                    for (i = k; i < m; i++) 
+                    {
+                        QR_(i,j) += s*QR_(i,k);
+                    }
+                }
+            }
+            Rdiag(k) = -nrm;
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    long qr_decomposition<matrix_exp_type>::
+    nr (
+    ) const
+    {
+        return m;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    long qr_decomposition<matrix_exp_type>::
+    nc (
+    ) const
+    {
+        return n;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    bool qr_decomposition<matrix_exp_type>::
+    is_full_rank(
+    ) const		
+    {
+        type eps = max(abs(Rdiag));
+        if (eps != 0)
+            eps *= std::sqrt(std::numeric_limits<type>::epsilon())/100;
+        else
+            eps = 1;  // there is no max so just use 1
+
+        // check if any of the elements of Rdiag are effectively 0
+        return min(abs(Rdiag)) > eps;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename qr_decomposition<matrix_exp_type>::matrix_type qr_decomposition<matrix_exp_type>::
+    get_householder (
+    )  const
+    {
+        return lowerm(QR_);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename qr_decomposition<matrix_exp_type>::matrix_type qr_decomposition<matrix_exp_type>::
+    get_r(
+    ) const
+    {
+        matrix_type R(n,n);
+        for (long i = 0; i < n; i++) 
+        {
+            for (long j = 0; j < n; j++) 
+            {
+                if (i < j) 
+                {
+                    R(i,j) = QR_(i,j);
+                } 
+                else if (i == j) 
+                {
+                    R(i,j) = Rdiag(i);
+                } 
+                else 
+                {
+                    R(i,j) = 0.0;
+                }
+            }
+        }
+        return R;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    const typename qr_decomposition<matrix_exp_type>::matrix_type qr_decomposition<matrix_exp_type>::
+    get_q(
+    ) const
+    {
+        long i=0, j=0, k=0;
+
+        matrix_type Q(m,n);
+        for (k = n-1; k >= 0; k--) 
+        {
+            for (i = 0; i < m; i++) 
+            {
+                Q(i,k) = 0.0;
+            }
+            Q(k,k) = 1.0;
+            for (j = k; j < n; j++) 
+            {
+                if (QR_(k,k) != 0) 
+                {
+                    type s = 0.0;
+                    for (i = k; i < m; i++) 
+                    {
+                        s += QR_(i,k)*Q(i,j);
+                    }
+                    s = -s/QR_(k,k);
+                    for (i = k; i < m; i++) 
+                    {
+                        Q(i,j) += s*QR_(i,k);
+                    }
+                }
+            }
+        }
+        return Q;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    template <typename EXP>
+    const typename qr_decomposition<matrix_exp_type>::matrix_type qr_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(B.nr() == nr(),
+            "\tconst matrix_type qr_decomposition::solve(B)"
+            << "\n\tInvalid inputs were given to this function"
+            << "\n\tB.nr():         " << B.nr()
+            << "\n\tnr():           " << nr()
+            << "\n\tthis:           " << this
+            );
+
+        // just call the right version of the solve function
+        if (B.nc() == 1)
+            return solve_vect(B);
+        else
+            return solve_mat(B);
+    }
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+//                           Private member functions
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    template <typename EXP>
+    const typename qr_decomposition<matrix_exp_type>::matrix_type qr_decomposition<matrix_exp_type>::
+    solve_vect(
+        const matrix_exp<EXP>& B
+    ) const
+    {
+
+        column_vector_type x(B);
+
+        // Compute Y = transpose(Q)*B
+        for (long k = 0; k < n; k++) 
+        {
+            type s = 0.0; 
+            for (long i = k; i < m; i++) 
+            {
+                s += QR_(i,k)*x(i);
+            }
+            s = -s/QR_(k,k);
+            for (long i = k; i < m; i++) 
+            {
+                x(i) += s*QR_(i,k);
+            }
+        }
+        // Solve R*X = Y;
+        for (long k = n-1; k >= 0; k--) 
+        {
+            x(k) /= Rdiag(k);
+            for (long i = 0; i < k; i++) 
+            {
+                x(i) -= x(k)*QR_(i,k);
+            }
+        }
+
+
+        /* return n x 1 portion of x */
+        return colm(x,0,n);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_exp_type>
+    template <typename EXP>
+    const typename qr_decomposition<matrix_exp_type>::matrix_type qr_decomposition<matrix_exp_type>::
+    solve_mat(
+        const matrix_exp<EXP>& B
+    ) const
+    {
+        const long nx = B.nc(); 
+        matrix_type X(B);
+        long i=0, j=0, k=0;
+
+        // Compute Y = transpose(Q)*B
+        for (k = 0; k < n; k++) 
+        {
+            for (j = 0; j < nx; j++) 
+            {
+                type s = 0.0; 
+                for (i = k; i < m; i++) 
+                {
+                    s += QR_(i,k)*X(i,j);
+                }
+                s = -s/QR_(k,k);
+                for (i = k; i < m; i++) 
+                {
+                    X(i,j) += s*QR_(i,k);
+                }
+            }
+        }
+        // Solve R*X = Y;
+        for (k = n-1; k >= 0; k--) 
+        {
+            for (j = 0; j < nx; j++) 
+            {
+                X(k,j) /= Rdiag(k);
+            }
+            for (i = 0; i < k; i++) 
+            {
+                for (j = 0; j < nx; j++) 
+                {
+                    X(i,j) -= X(k,j)*QR_(i,k);
+                }
+            }
+        }
+
+        /* return n x nx portion of X */
+        return subm(X,0,0,n,nx);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+
+} 
+
+#endif // DLIB_MATRIX_QR_DECOMPOSITION_H 
+
+
+

dlib/matrix/matrix_utilities.h

@@ -2030,13 +2030,13 @@ convergence:
     template <
         typename EXP
         >
-    inline const typename matrix_exp<EXP>::matrix_type cholesky_decomposition (
+    inline const typename matrix_exp<EXP>::matrix_type chol (
         const matrix_exp<EXP>& A
     )
     {
         DLIB_ASSERT(A.nr() == A.nc(), 
-            "\tconst matrix cholesky_decomposition(const matrix_exp& A)"
-            << "\n\tYou can only apply the cholesky_decomposition to a square matrix"
+            "\tconst matrix chol(const matrix_exp& A)"
+            << "\n\tYou can only apply the chol to a square matrix"
             << "\n\tA.nr(): " << A.nr()
             << "\n\tA.nc(): " << A.nc() 
             );

dlib/matrix/matrix_utilities_abstract.h

@@ -725,7 +725,7 @@ namespace dlib
 
 // ----------------------------------------------------------------------------------------
 
-    const matrix_exp::matrix_type cholesky_decomposition (
+    const matrix_exp::matrix_type chol (
         const matrix_exp& A
     );
     /*!

dlib/test/CMakeLists.txt

@@ -44,6 +44,7 @@ set (tests
    matrix.cpp
    matrix2.cpp
    matrix3.cpp
+   matrix_la.cpp
    md5.cpp
    member_function_pointer.cpp
    metaprogramming.cpp

dlib/test/matrix.cpp

@@ -608,7 +608,7 @@ namespace
             // make m be symmetric
             m = m*trans(m);
 
-            matrix<double> L = cholesky_decomposition(m);
+            matrix<double> L = chol(m);
             DLIB_CASSERT(equal(L*trans(L), m), "");
 
             DLIB_CASSERT(equal(inv(m), inv_upper_triangular(trans(L))*inv_lower_triangular((L))), "") 
@@ -633,7 +633,7 @@ namespace
             // make m be symmetric
             m = m*trans(m);
 
-            matrix<double> L = cholesky_decomposition(m);
+            matrix<double> L = chol(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))), "") 
@@ -667,7 +667,7 @@ namespace
 
             matrix<double,6,6> m(identity_matrix<double>(6)*4.5);
 
-            matrix<double> L = cholesky_decomposition(m);
+            matrix<double> L = chol(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/test/matrix_la.cpp

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+
+
+#include <dlib/matrix.h>
+#include <sstream>
+#include <string>
+#include <cstdlib>
+#include <ctime>
+#include <vector>
+#include "../stl_checked.h"
+#include "../array.h"
+#include "../rand.h"
+#include <dlib/string.h>
+
+#include "tester.h"
+
+namespace  
+{
+
+    using namespace test;
+    using namespace dlib;
+    using namespace std;
+
+    logger dlog("test.matrix_la");
+
+    dlib::rand::float_1a rnd;
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_type>
+    const typename matrix_type::matrix_type symm(const matrix_type& m) { return m*trans(m); }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename type>
+    const matrix<type> randm(long r, long c)
+    {
+        matrix<type> m(r,c);
+        for (long row = 0; row < m.nr(); ++row)
+        {
+            for (long col = 0; col < m.nc(); ++col)
+            {
+                m(row,col) = static_cast<type>(rnd.get_random_double()); 
+            }
+        }
+
+        return m;
+    }
+
+    template <typename type, long NR, long NC>
+    const matrix<type,NR,NC> randm()
+    {
+        matrix<type,NR,NC> m;
+        for (long row = 0; row < m.nr(); ++row)
+        {
+            for (long col = 0; col < m.nc(); ++col)
+            {
+                m(row,col) = static_cast<type>(rnd.get_random_double()); 
+            }
+        }
+
+        return m;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_type>
+    void test_qr ( const matrix_type& m)
+    {
+        typedef typename matrix_type::type type;
+        const type eps = 10*max(abs(m))*sqrt(std::numeric_limits<type>::epsilon());
+        dlog << LDEBUG << "test_qr():  " << m.nr() << " x " << m.nc() << "  eps: " << eps;
+        print_spinner();
+
+
+        qr_decomposition<matrix_type> test(m);
+
+
+        DLIB_CASSERT(test.nr() == m.nr(),"");
+        DLIB_CASSERT(test.nc() == m.nc(),"");
+
+
+        type temp;
+        DLIB_CASSERT( (temp= max(abs(test.get_q()*test.get_r() - m))) < eps,temp);
+
+        // none of the matrices we should be passing in to test_qr() should be non-full rank.  
+        DLIB_CASSERT(test.is_full_rank() == true,"");
+
+        if (m.nr() == m.nc())
+        {
+            matrix<type> m2;
+            matrix<type,0,1> col;
+
+            m2 = identity_matrix<type>(m.nr());
+            DLIB_CASSERT(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2)));
+            m2 = randm<type>(m.nr(),5);
+            DLIB_CASSERT(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2)));
+            m2 = randm<type>(m.nr(),1);
+            DLIB_CASSERT(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2)));
+            col = randm<type>(m.nr(),1);
+            DLIB_CASSERT(equal(m*test.solve(col), col,eps),max(abs(m*test.solve(m2)- m2)));
+        }
+        else
+        {
+            DLIB_CASSERT(dlib::equal(pinv(m), test.solve(identity_matrix<type>(m.nr())), eps), 
+                        max(abs(pinv(m) - test.solve(identity_matrix<type>(m.nr())))) );
+        }
+
+        // now make us a non-full rank matrix
+        if (m.nc() > 1)
+        {
+            matrix<type> sm(m);
+            set_colm(sm,0) = colm(sm,1);
+
+            qr_decomposition<matrix_type> test2(sm);
+            DLIB_CASSERT( (temp= max(abs(test.get_q()*test.get_r() - m))) < eps,temp);
+
+            if (test2.nc() < 100)
+            {
+                DLIB_CASSERT(test2.is_full_rank() == false,"eps: " << eps);
+            }
+
+        }
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_type>
+    void test_lu ( const matrix_type& m)
+    {
+        typedef typename matrix_type::type type;
+        const type eps = 10*max(abs(m))*sqrt(std::numeric_limits<type>::epsilon());
+        dlog << LDEBUG << "test_lu():  " << m.nr() << " x " << m.nc() << "  eps: " << eps;
+        print_spinner();
+
+
+        lu_decomposition<matrix_type> test(m);
+
+        DLIB_CASSERT(test.is_square() == (m.nr() == m.nc()), "");
+
+        DLIB_CASSERT(test.nr() == m.nr(),"");
+        DLIB_CASSERT(test.nc() == m.nc(),"");
+
+        type temp;
+        DLIB_CASSERT( (temp= max(abs(test.get_l()*test.get_u() - rowm(m,test.get_pivot())))) < eps,temp);
+
+        if (test.is_square())
+        {
+            // none of the matrices we should be passing in to test_lu() should be singular.  
+            DLIB_CASSERT (abs(test.det()) > eps/100, "det: " << test.det() );
+            dlog << LDEBUG << "big det: " << test.det();
+
+            DLIB_CASSERT(test.is_singular() == false,"");
+
+            matrix<type> m2;
+            matrix<type,0,1> col;
+
+            m2 = identity_matrix<type>(m.nr());
+            DLIB_CASSERT(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2)));
+            m2 = randm<type>(m.nr(),5);
+            DLIB_CASSERT(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2)));
+            m2 = randm<type>(m.nr(),1);
+            DLIB_CASSERT(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2)));
+            col = randm<type>(m.nr(),1);
+            DLIB_CASSERT(equal(m*test.solve(col), col,eps),max(abs(m*test.solve(m2)- m2)));
+
+            // now make us a singular matrix
+            if (m.nr() > 1)
+            {
+                matrix<type> sm(m);
+                set_colm(sm,0) = colm(sm,1);
+
+                lu_decomposition<matrix_type> test2(sm);
+                DLIB_CASSERT( (temp= max(abs(test2.get_l()*test2.get_u() - rowm(sm,test2.get_pivot())))) < eps,temp);
+
+                // these checks are only accurate for small matrices
+                if (test2.nr() < 100)
+                {
+                    DLIB_CASSERT(test2.is_singular() == true,"det: " << test2.det());
+                    DLIB_CASSERT(abs(test2.det()) < eps,"det: " << test2.det());
+                }
+
+            }
+        }
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_type>
+    void test_cholesky ( const matrix_type& m)
+    {
+        typedef typename matrix_type::type type;
+        const type eps = 10*max(abs(m))*sqrt(std::numeric_limits<type>::epsilon());
+        dlog << LDEBUG << "test_cholesky():  " << m.nr() << " x " << m.nc() << "  eps: " << eps;
+        print_spinner();
+
+
+        cholesky_decomposition<matrix_type> test(m);
+
+        // none of the matrices we should be passing in to test_cholesky() should be non-spd.  
+        DLIB_CASSERT(test.is_spd() == true,  "");
+
+        type temp;
+        DLIB_CASSERT( (temp= max(abs(test.get_l()*trans(test.get_l()) - m))) < eps,temp);
+
+
+        matrix<type> m2;
+        matrix<type,0,1> col;
+
+        m2 = identity_matrix<type>(m.nr());
+        DLIB_CASSERT(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2)));
+        m2 = randm<type>(m.nr(),5);
+        DLIB_CASSERT(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2)));
+        m2 = randm<type>(m.nr(),1);
+        DLIB_CASSERT(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2)));
+        col = randm<type>(m.nr(),1);
+        DLIB_CASSERT(equal(m*test.solve(col), col,eps),max(abs(m*test.solve(m2)- m2)));
+
+        // now make us a non-spd matrix
+        if (m.nr() > 1)
+        {
+            matrix<type> sm(lowerm(m));
+            sm(1,1) = 0;
+
+            cholesky_decomposition<matrix_type> test2(sm);
+            DLIB_CASSERT(test2.is_spd() == false,  test2.get_l());
+
+
+            cholesky_decomposition<matrix_type> test3(sm*trans(sm));
+            DLIB_CASSERT(test3.is_spd() == false,  test3.get_l());
+
+            sm = sm*trans(sm);
+            sm(1,1) = 5;
+            sm(1,0) -= 1;
+            cholesky_decomposition<matrix_type> test4(sm);
+            DLIB_CASSERT(test4.is_spd() == false,  test4.get_l());
+        }
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_type>
+    void test_eigenvalue ( const matrix_type& m)
+    {
+        typedef typename matrix_type::type type;
+        const type eps = max(abs(m))*sqrt(std::numeric_limits<type>::epsilon());
+        dlog << LDEBUG << "test_eigenvalue():  " << m.nr() << " x " << m.nc() << "  eps: " << eps;
+        print_spinner();
+
+
+        eigenvalue_decomposition<matrix_type> test(m);
+
+        DLIB_CASSERT(test.dim() == m.nr(), "");
+
+        // make sure all the various ways of asking for the eigenvalues are actually returning a
+        // consistent set of eigenvalues.
+        DLIB_CASSERT(equal(real(test.get_eigenvalues()), test.get_real_eigenvalues(), eps), ""); 
+        DLIB_CASSERT(equal(imag(test.get_eigenvalues()), test.get_imag_eigenvalues(), eps), ""); 
+        DLIB_CASSERT(equal(real(diag(test.get_d())), test.get_real_eigenvalues(), eps), ""); 
+        DLIB_CASSERT(equal(imag(diag(test.get_d())), test.get_imag_eigenvalues(), eps), ""); 
+
+        const matrix<type> V = test.get_pseudo_v();
+        const matrix<type> D = test.get_pseudo_d();
+        const matrix<complex<type> > CV = test.get_v();
+        const matrix<complex<type> > CD = test.get_d();
+        const matrix<complex<type> > CM = complex_matrix(m, uniform_matrix<type>(m.nr(),m.nc(),0));
+
+        DLIB_CASSERT(V.nr() == test.dim(),"");
+        DLIB_CASSERT(V.nc() == test.dim(),"");
+        DLIB_CASSERT(D.nr() == test.dim(),"");
+        DLIB_CASSERT(D.nc() == test.dim(),"");
+
+        // CD is a diagonal matrix
+        DLIB_CASSERT(diagm(diag(CD)) == CD,"");
+
+        // verify that these things are actually eigenvalues and eigenvectors of m
+        DLIB_CASSERT(max(abs(m*V - V*D)) < eps, "");
+        DLIB_CASSERT(max(norm(CM*CV - CV*CD)) < eps, "");
+
+        // if m is a symmetric matrix
+        if (max(abs(m-trans(m))) < 1e-5)
+        {
+            dlog << LTRACE << "m is symmetric";
+            // there aren't any imaginary eigenvalues 
+            DLIB_CASSERT(max(abs(test.get_imag_eigenvalues())) < eps, ""); 
+            DLIB_CASSERT(diagm(diag(D)) == D,"");
+
+            // V is orthogonal
+            DLIB_CASSERT(equal(V*trans(V), identity_matrix<type>(test.dim()), eps), "");
+            DLIB_CASSERT(equal(m , V*D*trans(V), eps), "");
+        }
+        else
+        {
+            dlog << LTRACE << "m is NOT symmetric";
+            DLIB_CASSERT(equal(m , V*D*inv(V), eps), max(abs(m - V*D*inv(V))));
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    void matrix_test_double()
+    {
+
+    // -------------------------------
+
+        test_lu(10*randm<double>(1,1));
+        test_lu(10*randm<double>(2,2));
+        test_lu(10*symm(randm<double>(2,2)));
+        test_lu(10*randm<double>(4,4));
+        test_lu(10*randm<double>(9,4));
+        test_lu(10*randm<double>(3,8));
+        test_lu(10*randm<double>(15,15));
+        test_lu(2*symm(randm<double>(15,15)));
+        test_lu(10*randm<double>(100,100));
+        test_lu(10*randm<double>(137,200));
+        test_lu(10*randm<double>(200,101));
+
+        test_lu(10*randm<double,1,1>());
+        test_lu(10*randm<double,2,2>());
+        test_lu(10*randm<double,4,3>());
+        test_lu(10*randm<double,4,4>());
+        test_lu(10*randm<double,9,4>());
+        test_lu(10*randm<double,3,8>());
+        test_lu(10*randm<double,15,15>());
+        test_lu(10*randm<double,100,100>());
+        test_lu(10*randm<double,137,200>());
+        test_lu(10*randm<double,200,101>());
+
+    // -------------------------------
+
+        test_cholesky(uniform_matrix<double>(1,1,1) + 10*symm(randm<double>(1,1)));
+        test_cholesky(uniform_matrix<double>(2,2,1) + 10*symm(randm<double>(2,2)));
+        test_cholesky(uniform_matrix<double>(3,3,1) + 10*symm(randm<double>(3,3)));
+        test_cholesky(uniform_matrix<double>(4,4,1) + 10*symm(randm<double>(4,4)));
+        test_cholesky(uniform_matrix<double>(15,15,1) + 10*symm(randm<double>(15,15)));
+        test_cholesky(uniform_matrix<double>(101,101,1) + 10*symm(randm<double>(101,101)));
+
+    // -------------------------------
+
+        test_qr(10*randm<double>(1,1));
+        test_qr(10*randm<double>(2,2));
+        test_qr(10*symm(randm<double>(2,2)));
+        test_qr(10*randm<double>(4,4));
+        test_qr(10*randm<double>(9,4));
+        test_qr(10*randm<double>(15,15));
+        test_qr(2*symm(randm<double>(15,15)));
+        test_qr(10*randm<double>(100,100));
+        test_qr(10*randm<double>(237,200));
+        test_qr(10*randm<double>(200,101));
+
+        test_qr(10*randm<double,1,1>());
+        test_qr(10*randm<double,2,2>());
+        test_qr(10*randm<double,4,3>());
+        test_qr(10*randm<double,4,4>());
+        test_qr(10*randm<double,9,4>());
+        test_qr(10*randm<double,15,15>());
+        test_qr(10*randm<double,100,100>());
+
+    // -------------------------------
+
+        test_eigenvalue(10*randm<double>(1,1));
+        test_eigenvalue(10*randm<double>(2,2));
+        test_eigenvalue(10*randm<double>(3,3));
+        test_eigenvalue(10*randm<double>(4,4));
+        test_eigenvalue(10*randm<double>(15,15));
+        test_eigenvalue(10*randm<double>(150,150));
+
+        test_eigenvalue(10*randm<double,1,1>());
+        test_eigenvalue(10*randm<double,2,2>());
+        test_eigenvalue(10*randm<double,3,3>());
+
+        test_eigenvalue(10*symm(randm<double>(1,1)));
+        test_eigenvalue(10*symm(randm<double>(2,2)));
+        test_eigenvalue(10*symm(randm<double>(3,3)));
+        test_eigenvalue(10*symm(randm<double>(4,4)));
+        test_eigenvalue(10*symm(randm<double>(15,15)));
+        test_eigenvalue(10*symm(randm<double>(150,150)));
+
+    // -------------------------------
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    void matrix_test_float()
+    {
+
+    // -------------------------------
+
+        test_lu(3*randm<float>(1,1));
+        test_lu(3*randm<float>(2,2));
+        test_lu(3*randm<float>(4,4));
+        test_lu(3*randm<float>(9,4));
+        test_lu(3*randm<float>(3,8));
+        test_lu(3*randm<float>(137,200));
+        test_lu(3*randm<float>(200,101));
+
+        test_lu(3*randm<float,1,1>());
+        test_lu(3*randm<float,2,2>());
+        test_lu(3*randm<float,4,3>());
+        test_lu(3*randm<float,4,4>());
+        test_lu(3*randm<float,9,4>());
+        test_lu(3*randm<float,3,8>());
+        test_lu(3*randm<float,137,200>());
+        test_lu(3*randm<float,200,101>());
+
+    // -------------------------------
+
+        test_cholesky(uniform_matrix<float>(1,1,1) + 2*symm(randm<float>(1,1)));
+        test_cholesky(uniform_matrix<float>(2,2,1) + 2*symm(randm<float>(2,2)));
+        test_cholesky(uniform_matrix<float>(3,3,1) + 2*symm(randm<float>(3,3)));
+
+    // -------------------------------
+
+        test_qr(3*randm<float>(1,1));
+        test_qr(3*randm<float>(2,2));
+        test_qr(3*randm<float>(4,4));
+        test_qr(3*randm<float>(9,4));
+        test_qr(3*randm<float>(237,200));
+
+        test_qr(3*randm<float,1,1>());
+        test_qr(3*randm<float,2,2>());
+        test_qr(3*randm<float,4,3>());
+        test_qr(3*randm<float,4,4>());
+        test_qr(3*randm<float,9,4>());
+
+    // -------------------------------
+
+
+        test_eigenvalue(10*randm<float>(1,1));
+        test_eigenvalue(10*randm<float>(2,2));
+        test_eigenvalue(10*randm<float>(3,3));
+        test_eigenvalue(10*randm<float>(4,4));
+        test_eigenvalue(10*randm<float>(15,15));
+        test_eigenvalue(10*randm<float>(150,150));
+
+        test_eigenvalue(10*randm<float,1,1>());
+        test_eigenvalue(10*randm<float,2,2>());
+        test_eigenvalue(10*randm<float,3,3>());
+
+        test_eigenvalue(10*symm(randm<float>(1,1)));
+        test_eigenvalue(10*symm(randm<float>(2,2)));
+        test_eigenvalue(10*symm(randm<float>(3,3)));
+        test_eigenvalue(10*symm(randm<float>(4,4)));
+        test_eigenvalue(10*symm(randm<float>(15,15)));
+        test_eigenvalue(10*symm(randm<float>(150,150)));
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    class matrix_tester : public tester
+    {
+    public:
+        matrix_tester (
+        ) :
+            tester ("test_matrix_la",
+                    "Runs tests on the matrix component.")
+        {
+            rnd.set_seed(cast_to_string(time(0)));
+        }
+
+        void perform_test (
+        )
+        {
+            dlog << LINFO << "seed string: " << rnd.get_seed();
+
+            dlog << LINFO << "begin testing with double";
+            matrix_test_double();
+            dlog << LINFO << "begin testing with float";
+            matrix_test_float();
+        }
+    } a;
+
+}
+
+
+