Skip to content

Commit

Permalink
Implemented determinant() + a bunch of other stuff in Matrix
Browse files Browse the repository at this point in the history
  • Loading branch information
@hosseinmoein
hosseinmoein committed Jan 9, 2025
1 parent 76e93d1 commit 6705f01
Show file tree
Hide file tree
Showing 3 changed files with 350 additions and 5 deletions.
89 changes: 89 additions & 0 deletions include/DataFrame/Utils/Matrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,6 +281,95 @@ class Matrix {
template<typename MA1, typename MA2, typename MA3>
inline void svd(MA1 &U, MA2 &S, MA3 &V, bool full_size = true) const;

// In linear algebra, the LU decomposition is a matrix decomposition
// which writes a matrix as the product of a lower and upper triangular
// amtrices. This decomposition is used in numerical analysis to solve
// systems of linear equations or calculate the determinant.
//
// LU decomposition of a square matrix A is a decomposition of A as
// A = LU
//
// where L and U are lower and upper triangular matrices (of the same
// size), respectively. This means that L has only zeros above the
// diagonal and U has only zeros below the diagonal.
//
// The LU decomposition is twice as efficient as the QR decomposition.
//
// NOTE: The LU decompostion with pivoting always exists, even if the
// matrix is singular. 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 false.
//
template<typename MA1, typename MA2>
inline void lud(MA1 &L, MA2 &U) const;

// In 2-D, when you talk about the point (2, 4), you can think of the
// 2 and 4 as directions to get from the origin to the point
// (move 2 units in the x direction and 4 in the y direction). In
// a 3-D system, the same idea holds - (1, 3, 7) means start at the
// origin (0, 0, 0), go 1 unit in the x direction, 3 in the y direction,
// and 7 in the z direction.
// The determinant of a 1x1 matrix is the signed length of the line
// from the origin to the point. It's positive if the point is in the
// positive x direction, negative if in the other direction.
// In 2-D, look at the matrix as two 2-dimensional points on the plane,
// and complete the parallelogram that includes those two points and
// the origin. The (signed) area of this parallelogram is the
// determinant. If you sweep clockwise from the first to the second,
// the determinant is negative, otherwise, positive.
// In 3-D, look at the matrix as 3 3-dimensional points in space.
// Complete the parallepipe that includes these points and the origin,
// and the determinant is the (signed) volume of the parallelepipe.
// The same idea works in any number of dimensions. The determinant
// is just the (signed) volume of the n-dimensional parallelepipe.
// Note that length, area, volume are the _volumes_ in 1, 2, and
// 3-dimensional spaces. A similar concept of volume exists for
// Euclidean space of any dimensionality.
//
// [HM]: A matrix that has determinant equals zero is singular.
// It means it has at least 2 vectors that are linearly
// dependent, therefore they cannot enclose an area or a space.
//
// For how determinant of N X N matrices are calulated
// see: http://mathworld.wolfram.com/Determinant.html
// For Laplace Expension
// see: http://en.wikipedia.org/wiki/Cofactor_expansion
//
// NOTE: This is a relatively _expensive_ calculation. Its complexity
// is O(n!), where n is the number of rows or columns.
//
inline value_type determinant() const;

// Minor of a matrix is the same matrix with the specified row
// and column taken out.
//
// NOTE: Linear algebrically speaking, a minor of a matrix is the
// _detereminant_ of some smaller square matrix, cut down from
// the original matrix.
//
template<typename MA>
inline void
get_minor(MA &mmatrix, size_type drow, size_type dcol) const noexcept;

// A Cofactor of a matrix is the determinant of a minor of the matrix.
// You can also say cofactor is the signed minor of the matrix.
//
inline value_type cofactor(size_type row, size_type column) const;

// The Adjoint of a matrix is formed by taking the transpose of the
// cofactors matrix of the original matrix.
//
template<typename MA>
inline void adjoint(MA &amatrix) const;

// A "centered matrix" is a matrix where the mean of each column has been
// adjusted to be zero, meaning that the data within each column is
// centered around the value of zero
//
inline void center() noexcept;
template<typename MA>
inline void get_centered(MA &cmatrix) const noexcept;

private:

static constexpr size_type NOPOS_ = static_cast<size_type>(-9);
Expand Down
231 changes: 226 additions & 5 deletions include/DataFrame/Utils/Matrix.tcc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -320,10 +320,10 @@ Matrix<T, MO, IS_SYM>::transpose2() const noexcept {
// ----------------------------------------------------------------------------

template<typename T, matrix_orient MO, bool IS_SYM>
inline Matrix<T, MO, IS_SYM>::size_type
Matrix<T, MO, IS_SYM>::ppivot_(size_type pivot_row,
size_type self_rows,
size_type self_cols) noexcept {
inline Matrix<T, MO, IS_SYM>::size_type Matrix<T, MO, IS_SYM>::
ppivot_(size_type pivot_row,
size_type self_rows,
size_type self_cols) noexcept {

size_type max_row { pivot_row };
value_type max_value { value_type(std::fabs(at(pivot_row, pivot_row))) };
Expand Down Expand Up @@ -1332,7 +1332,6 @@ eigen_space(MA1 &eigenvalues, MA2 &eigenvectors, bool sort_values) const {
throw DataFrameError("Matrix::eigen_space(): Matrix must be squared");
#endif // HMDF_SANITY_EXCEPTIONS


MA1 tmp_evals { 1, cols() };
MA2 tmp_evecs { rows(), cols() };
Matrix imagi { 1, cols() }; // Imaginary part
Expand Down Expand Up @@ -1886,6 +1885,228 @@ svd(MA1 &U, MA2 &S, MA3 &V, bool full_size) const {
return;
}

// ----------------------------------------------------------------------------

template<typename T, matrix_orient MO, bool IS_SYM>
template<typename MA1, typename MA2>
void Matrix<T, MO, IS_SYM>::
lud(MA1 &L, MA2 &U) const {

#ifdef HMDF_SANITY_EXCEPTIONS
if (! is_square())
throw DataFrameError("Matrix::lud(): Matrix must be squared");
#endif // HMDF_SANITY_EXCEPTIONS

Matrix tmp = *this;

for (size_type c = 0; c < cols(); ++c) {
// Find pivot.
//
size_type p { c };

for (size_type r = c + 1; r < rows(); ++r)
if (tmp(r, c) < tmp(p, c))
p = r;

// Exchange if necessary.
//
if (p != c)
for (size_type cc = 0; cc < cols(); ++cc)
std::swap(tmp(p, cc), tmp(c, cc));

// Compute multipliers and eliminate c-th column.
//
if (tmp(c, c) != value_type(0))
for (size_type r = c + 1; r < rows(); ++r) {
tmp(r, c) /= tmp(c, c);

for (size_type cc = c + 1; cc < cols(); ++cc)
tmp(r, cc) -= tmp(r, c) * tmp(c, cc);
}
}

MA1 l_tmp { rows(), cols(), 0 };

for (size_type r = 0; r < rows(); ++r) {
for (size_type c = 0; c < cols(); ++c) {
if (r > c)
l_tmp(r, c) = tmp(r, c);
else if (r == c)
l_tmp(r, c) = value_type(1);
}
}

MA2 u_tmp { rows(), cols(), 0 };

for (size_type c = 0; c < cols(); ++c) {
for (size_type r = 0; r < rows(); ++r) {
if (c <= r)
u_tmp(c, r) = tmp(c, r);
}
}

L.swap(l_tmp);
U.swap(u_tmp);

return;
}

// ----------------------------------------------------------------------------

template<typename T, matrix_orient MO, bool IS_SYM>
Matrix<T, MO, IS_SYM>::value_type
Matrix<T, MO, IS_SYM>::determinant() const {

#ifdef HMDF_SANITY_EXCEPTIONS
if (! is_square())
throw DataFrameError("Matrix::determinant(): Matrix must be squared");
#endif // HMDF_SANITY_EXCEPTIONS

const size_type sz = rows();
value_type result { 1 };

if (sz == 2) {
result = at(0, 0) * at(1, 1) - at(0, 1) * at(1, 0);
}
else if (sz < 8) { // Cofacter way
Matrix tmp = *this;

for (size_type r = 0; r < sz; ++r) {
const size_type indx = tmp.ppivot_(r, sz, sz);

if (indx == NOPOS_)
return (value_type(0));

if (indx != 0)
result = -result;

const value_type diag = tmp (r, r);

result *= diag;
for (size_type rr = r + 1; rr < sz; ++rr) {
const value_type piv = tmp (rr, r) / diag;

for (size_type c = r + 1; c < sz; ++c)
tmp (rr, c) -= piv * tmp (r, c);
}
}
}
else { // LUD way
Matrix<T, matrix_orient::row_major> L { sz, sz, 0 };
Matrix<T, matrix_orient::column_major> U { sz, sz, 0 };

lud(L, U);
for (size_type r = 0; r < sz; ++r)
result *= U(r, r);
}

return (result);
}

// ----------------------------------------------------------------------------

template<typename T, matrix_orient MO, bool IS_SYM>
template<typename MA>
void Matrix<T, MO, IS_SYM>::
get_minor (MA &mmatrix, size_type drow, size_type dcol) const noexcept {

mmatrix.resize(rows() - 1, cols() - 1, 0);

size_type mrow = 0;

for (size_type r = 0; r < rows(); ++r) {
if (r != drow) {
size_type mcol = 0;

for (size_type c = 0; c < cols(); ++c) {
if (c != dcol)
mmatrix(mrow, mcol++) = at(r, c);
}

mrow += 1; // Increase minor row-count
}
}

return;
}

// ----------------------------------------------------------------------------

template<typename T, matrix_orient MO, bool IS_SYM>
inline typename Matrix<T, MO, IS_SYM>::value_type
Matrix<T, MO, IS_SYM>::cofactor (size_type row, size_type column) const {

#ifdef HMDF_SANITY_EXCEPTIONS
if (! is_square())
throw DataFrameError("Matrix::cofactor(): Matrix must be squared");
#endif // HMDF_SANITY_EXCEPTIONS

Matrix<T, matrix_orient::row_major, IS_SYM> tmp;

get_minor(tmp, row, column);
return ((row + column) % 2 ? -tmp.determinant() : tmp.determinant());
}

// ----------------------------------------------------------------------------

template<typename T, matrix_orient MO, bool IS_SYM>
template<typename MA>
inline void Matrix<T, MO, IS_SYM>::adjoint (MA &that) const {

#ifdef HMDF_SANITY_EXCEPTIONS
if (! is_square())
throw DataFrameError("Matrix::adjoint(): Matrix must be squared");
#endif // HMDF_SANITY_EXCEPTIONS

that.resize(rows(), cols());
for (size_type r = 0; r < rows(); ++r)
for (size_type c = 0; c < cols(); ++c)
that(c, r) = cofactor(r, c);

return;
}

// ----------------------------------------------------------------------------

template<typename T, matrix_orient MO, bool IS_SYM>
void Matrix<T, MO, IS_SYM>::center() noexcept {

for (size_type c = 0; c < cols(); ++c) {
value_type mean { 0 };

for (size_type r = 0; r < rows(); ++r)
mean += at(r, c);
mean /= value_type(rows());

for (size_type r = 0; r < rows(); ++r)
at(r, c) -= mean;
}

return;
}

// ----------------------------------------------------------------------------

template<typename T, matrix_orient MO, bool IS_SYM>
template<typename MA>
void Matrix<T, MO, IS_SYM>::get_centered(MA &cmatrix) const noexcept {

cmatrix.resize(rows(), cols(), 0);
if constexpr (MO == matrix_orient::column_major) {
for (size_type c = 0; c < cols(); ++c)
for (size_type r = c; r < rows(); ++r)
cmatrix(r, c) = at(r, c);
}
else {
for (size_type r = 0; r < rows(); ++r)
for (size_type c = r; c < cols(); ++c)
cmatrix(r, c) = at(r, c);
}
cmatrix.center();

return;
}

} // namespace hmdf

// ----------------------------------------------------------------------------
Expand Down
35 changes: 35 additions & 0 deletions test/matrix_tester.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -499,6 +499,41 @@ int main(int, char *[]) {
assert(laplacian_mat(4, 4) == 2.0);
}

// Test Determinant and get_centered
//
{
using row_dmat_t = Matrix<double, matrix_orient::row_major>;

row_dmat_t mat1 { 2, 2 };
long value { 0 };

for (long r = 0; r < mat1.rows(); ++r)
for (long c = 0; c < mat1.cols(); ++c)
mat1(r, c) = double(++value);
assert(mat1.determinant() == -2.0);

row_dmat_t mat2 { 5, 5 };
std::vector<double> row0 = { 3.2, 11, 4.4, 3.3, 11 };
std::vector<double> row1 = { 4.5, 10.8, 8.3, 3.4, 12 };
std::vector<double> row2 = { 4.8, 9.1, 7.1, 3.6, 10.8 };
std::vector<double> row3 = { 5.1, 2, 7.8, 4.1, 1.1 };
std::vector<double> row4 = { 5.5, 0.5, 1.1, 4.12, 8 };

mat2.set_row(row0.begin(), 0);
mat2.set_row(row1.begin(), 1);
mat2.set_row(row2.begin(), 2);
mat2.set_row(row3.begin(), 3);
mat2.set_row(row4.begin(), 4);
assert((std::fabs(mat2.determinant() - 288.973) < 0.001));

row_dmat_t centered;

mat2.get_centered(centered);
assert((std::fabs(centered(0, 0) - 2.56) < 0.01));
assert((std::fabs(centered(2, 0) - -0.64) < 0.01));
assert((std::fabs(centered(4, 0) - -0.64) < 0.01));
}

test_thread_pool();
return (0);
}
Expand Down

0 comments on commit 6705f01

Please sign in to comment.