mat.cpp

这是个人脸识别程序

			*pSmallestPivot = fabs(LU(i,i));
	}

if (0 != (iRet = gsl_linalg_LU_solve(LU.m, pPerm, &bVec.vector, xVec)))
	SysErr("LUSolve gsl_linalg_LU_solve returned %d", iRet);

Vec x(xVec->data, nrows);	// convert gsl_vector to Vec
gsl_vector_free(xVec);
*ppPerm = pPerm;

return x;
}

//-----------------------------------------------------------------------------
// Same as above but with no LU or Perm params
// pSmallestPivot can be NULL, as above

Vec SolveWithLU (const Mat &A, Vec &b, double *pSmallestPivot)
{
#if GSL_RANGE_CHECK
CheckIsVector(b, "SolveWithLU b");
if (A.nrows() != b.nelems())
	SysErr("LUSolve b %dx%d doesn't conform to A %dx%d",
		b.nrows(), b.ncols(), A.nrows(), A.ncols());
#endif

int	Sign, iRet, nrows = A.nrows();

gsl_permutation *pPerm = gsl_permutation_alloc(nrows);
gsl_vector 		*xVec  = gsl_vector_alloc(nrows);
gsl_vector_view  bVec  = gsl_vector_view_array(b.m->data, nrows);
gsl_matrix		*pLU   = gsl_matrix_alloc(nrows, A.ncols());

gsl_matrix_memcpy(pLU, A.m);	// use copy not view because decomp destroys A

if (0 != (iRet = gsl_linalg_LU_decomp(pLU, pPerm, &Sign)))
	SysErr("LUSolve gsl_linalg_LU_decomp returned %d", iRet);

if (pSmallestPivot)		// find and store the smallest pivot?
	{
	*pSmallestPivot = GSL_POSINF;
	for (int i = A.nrows()-1; i >= 0; i--)	// smallest last so reduce nbr assignments
		if (fabs(gsl_matrix_get(pLU, i, i)) < *pSmallestPivot)
			*pSmallestPivot = fabs(gsl_matrix_get(pLU, i, i));
	}

if (0 != (iRet = gsl_linalg_LU_solve(pLU, pPerm, &bVec.vector, xVec)))
	SysErr("LUSolve gsl_linalg_LU_solve returned %d", iRet);

Vec x(xVec->data, nrows);	// convert gsl_vector to Vec
gsl_permutation_free(pPerm);
gsl_vector_free(xVec);
gsl_matrix_free(pLU);

return x;
}

//-----------------------------------------------------------------------------
// Like LUsolve but A must already be in LU form and pPerm must be initialized
// by a previous call to LUsolve.
//
// Calling this instead of LUSolve again is faster.

Vec SolveWithLUAgain (Mat &LU, Vec &b, gsl_permutation *pPerm)
{
#if GSL_RANGE_CHECK

CheckIsVector(b, "SolveWithLUAgain b");

if (LU.nrows() != b.nelems())
	SysErr("LUSolveAgain b %dx%d doesn't conform to LU %dx%d",
		b.nrows(), b.ncols(), LU.nrows(), LU.ncols());
#endif

int	iRet, nrows = LU.nrows();

gsl_vector *xVec = gsl_vector_alloc(nrows);

gsl_vector_view  bVec = gsl_vector_view_array(b.m->data, nrows);

if (0 != (iRet = gsl_linalg_LU_solve(LU.m, pPerm, &bVec.vector, xVec)))
	SysErr("LUSolveAgain gsl_linalg_LU_solve returned %d", iRet);

Vec x(xVec->data, nrows);	// convert gsl_vector to Vec
gsl_vector_free(xVec);

return x;
}

//-----------------------------------------------------------------------------
// A is not changed

Vec RefineLU (Mat &A, Mat &LU, gsl_permutation *pPerm, Vec &x, Vec &b)
{
int	iRet, nrows = A.nrows();

#if GSL_RANGE_CHECK
if (!b.fVector())
	SysErr("LURefine b %dx%d is not a vector", b.nrows(), b.ncols());

if (LU.nrows() != nrows || LU.ncols() != A.ncols())
	SysErr("LURefine A %dx%d EA %dx%d incompatible",
		nrows, A.ncols(), LU.nrows(), LU.ncols());

if (b.nrows() != nrows && b.ncols() != nrows)
	SysErr("LURefine neither b size %d %d != A.nrows %d", b.nrows(), b.ncols(), nrows);

if (pPerm->size != nrows)
	SysErr("LURefine pPerm->size %d != A.nrows %d", pPerm->size, nrows);
#endif

gsl_vector_view bVec = gsl_vector_view_array(b.m->data, nrows);
gsl_vector_view xVec = gsl_vector_view_array(x.m->data, nrows);
gsl_vector *resVec = gsl_vector_alloc(nrows);

if (0 != (iRet = gsl_linalg_LU_refine(A.m, LU.m, pPerm, &bVec.vector, &xVec.vector, resVec)))
	SysErr("LURefine gsl_linalg_LU_refine returned %d", iRet);

gsl_vector_free(resVec);

return x;
}

//-----------------------------------------------------------------------------
// This destroys A, or rather puts U into A

void SingularValDecomp (Mat &A, 		// io
			Mat &VTranspose, Vec &S)	// out
{
int iRet;

if (A.M() < A.N())	// check against GSL library limitation
	SysErr("SingularValDecomp: m < n not yet implemented A %dx%d", A.M(), A.N());

// if (A.N() == 1)		// check against GSL 1.4 library limitation, fixed in GSL 1.6
//	SysErr("SingularValDecomp: n == 1 not yet implemented A %dx%d", A.M(), A.N());

gsl_vector_view DVec = gsl_vector_view_array(S.m->data, A.N());
gsl_vector *W = gsl_vector_alloc(A.N());			// working vector

if (A.M() < (5 * A.N()) / 3)	// Golub & vanLoan p253 tells us the fastest SVD func to use
	{
	if (0 != (iRet = gsl_linalg_SV_decomp(A.m, VTranspose.m, &DVec.vector, W)))
		SysErr("SingularValDecomp gsl_linalg_SV_decomp returned %d", iRet);
	}
else
	{
	gsl_matrix *X = gsl_matrix_alloc(A.N(), A.N());		// working matrix

	if (0 != (iRet = gsl_linalg_SV_decomp_mod(A.m, X, VTranspose.m, &DVec.vector, W)))
		SysErr("SingularValDecomp gsl_linalg_SV_decomp_mod returned %d", iRet);

	gsl_matrix_free(X);
	}
gsl_vector_free(W);
}

//-----------------------------------------------------------------------------
// Solves Ax=b for x and returns x.
// A gets destroyed (actually it gets changed to the U in UDV).
//
// This returns the best answer in a least squares sense if no exact solution exists
// i.e if b is not in the column space of A or A is not invertible.

Vec SolveWithSVD (Mat &A, const Vec &b, bool fVerbose)
{
if (fVerbose)
	lprintf("SolveWithSVD %dx%d        ", A.nrows(), A.ncols());

int iRet;
Mat VTrans(A.N(), A.N());
Vec S(A.N());
Vec x(A.M());

SingularValDecomp(A, VTrans, S);

gsl_vector_view SVec = gsl_vector_view_array(S.m->data, A.N());
gsl_vector_view bVec = gsl_vector_view_array(b.m->data, A.M());
gsl_vector_view xVec = gsl_vector_view_array(x.m->data, A.M());

iRet = gsl_linalg_SV_solve(A.m, VTrans.m, &SVec.vector, &bVec.vector, &xVec.vector);

DASSERT(iRet == 0);

return x;
}

//-----------------------------------------------------------------------------
// returns true if the diag elements of square matrix A are strictly positive

static bool fStrictlyPositiveDiag (Mat &A)
{
for (int i = 0; i < A.nrows(); i++)
	if (A(i,i) <= 0)
		return false;
return true;
}

//-----------------------------------------------------------------------------
// Solves Ax=b for x and returns x.  Destroys A.
// This routine silently sets pfPosDef to false and immediately returns
// if A is not symmetric positive definite and it can't be solved by Cholesky.

Vec SolveWithCholesky (Mat &A, const Vec &b, bool fVerbose, bool *pfPosDef)
{
Vec x(A.nrows());

*pfPosDef = false;	// assume

if (0 == gsl_linalg_cholesky_decomp_return_on_EDOM(A.m))
	{
	int iRet;

	if (fVerbose)
		lprintf("SolveWithCholesky %dx%d   ", A.nrows(), A.ncols());

	gsl_vector_view bVec = gsl_vector_view_array(b.m->data, A.M());
	gsl_vector_view xVec = gsl_vector_view_array(x.m->data, A.M());

	if (0 != (iRet = gsl_linalg_cholesky_solve(A.m, &bVec.vector, &xVec.vector)))
		SysErr("SolveWithCholesky gsl_linalg_cholesky_solve returned %d", iRet);

	*pfPosDef = true;
	}
return x;
}

//-----------------------------------------------------------------------------
// Return true if A is positive definite

bool fPosDef (const Mat &A)
{
Mat ACopy(A); 	// SolveWithCholesky destroys A, so make a copy
Vec b(A.nrows());
bool fPosDef;
SolveWithCholesky(ACopy, b, false, &fPosDef);
return fPosDef;
}

//-----------------------------------------------------------------------------
// Like fPosDef() but uses the eig val method, and allows up to nAllowedZeroEigs eigs to near 0.
// Much slower than fPosDef() because has to find the eigs of A.

bool fPosDef1 (const Mat &A, int nAllowedZeroEigs)
{
int nelems = A.nrows();
Vec EigVals(nelems);
double Min = EigVals(0) / 1e6;

GetEigsForSymMat(A, EigVals, true, true);

for (int i = 0; i < nelems - nAllowedZeroEigs; i++)
	if (EigVals(i) <= Min)
		return false;

return true;
}

//-----------------------------------------------------------------------------
// Returns a matrix of eigenvectors for the symmetric matrix A.
// Eigenvectors not normalized. Also returns the eigen vals in EigenVals.
// If fSort is set, then eigenvecs and eigenvals are sorted, largest eigenval first.
// A must be a symmetric matrix. If fTestForSymmetry then SysErr if A not symmetric.
// A is not modified

Mat GetEigsForSymMat (const Mat &A, Vec &EigVals, bool fSort, bool fTestForSymmetry)
{
int iRet;

// some preliminary checks, GSL funcs will do more

DASSERT(A.nrows() == A.ncols());
DASSERT(EigVals.nrows() == A.nrows());

if (fTestForSymmetry && !fIsMatrixSymmetric(A, A(0,0) / 1e5))		//TODO magic 1e5 chosen based for ASM model building, make it a parameter?
	{
	//A.print("\nNon symmetric matrix\n", "%16.6g ");
	SysErr("GetEigsForSymMat matrix %dx%d not symmetric", A.nrows(), A.ncols());
	}
Mat EigVecs(A.ncols(), A.ncols());

Mat Work(A);		// gsl_eigen_symmv() destoys the A->m.data, so we use a temporary variable

gsl_eigen_symmv_workspace *W = gsl_eigen_symmv_alloc(Work.ncols());
gsl_vector_view EigValsVec = gsl_vector_view_array(EigVals.m->data, EigVals.nrows());

if (0 != (iRet = gsl_eigen_symmv(Work.m, &EigValsVec.vector, EigVecs.m, W)))
	SysErr("GetEigsForSymMat matrix %dx%d gsl_eigen_symmv returned %d", Work.nrows(), Work.ncols(), iRet);

gsl_eigen_symmv_free(W);

// if asked to, sort the eigenvalues and corresponding eigenvectors

if (fSort && 0 != (iRet =
		gsl_eigen_symmv_sort(&EigValsVec.vector, EigVecs.m, GSL_EIGEN_SORT_VAL_DESC)))
	{
	SysErr("GetEigsForSymMat matrix %dx%d gsl_eigen_symmv_sort returned %d", Work.nrows(), Work.ncols(), iRet);
	}
return EigVecs;
}

//-----------------------------------------------------------------------------
Mat NormalizeCols (const Mat &A)	// returns the mat but with cols normalized (length 1)
{
for (int iCol = 0; iCol < A.ncols(); iCol++)
	{
	MatView col = A.col(iCol);
	double norm = col.len();
	col /= norm;
	}
return A;
}

//-----------------------------------------------------------------------------
Mat HilbertMat (size_t n)		// returns a hilbert mat of the given size
{
Mat result(n, n);

for (int iRow = 0; iRow < n; iRow++)
	for (int iCol = 0; iCol < n; iCol++)
		gsl_matrix_set(result.m, iRow, iCol, 1.0 / (iRow + iCol + 1.0));

return result;
}

//-----------------------------------------------------------------------------
Mat IdentMat (size_t n)		// returns an identity mat of the given size
{
Mat result(n, n);
result.diagMe(1);
return result;
}

//-----------------------------------------------------------------------------
Vec AllOnesMat (size_t n, const char *sIsRowVec)	// returns an all ones mat of the given size
{
Vec result(n, sIsRowVec);
result.fill(1);
return result;
}

//-----------------------------------------------------------------------------
// This returns x * A * x.
// x is a vector.
// A is assumed to be a symmetric matrix (but only the upper right triangle of A is actually used).
//
// This function is equivalent to OneElemToDouble(x * A * x.t()), but is
// optimized for speed and is much faster.
//
// It's faster because we use the fact that A is symmetric to roughly
// halve the number of operations.
// We also bypass the overhead of normal matrix routines, which can be slow.

double xAx (const Vec &x, const Mat &A)
{
double xTemp[MAX_TEMP_ARRAY_LEN];	// copy x into here because it's faster accessing an array than a Vec

const int nelems = x.nelems();

ASSERT(nelems < MAX_TEMP_ARRAY_LEN);
DASSERT(A.nrows() == nelems);
DASSERT(A.ncols() == nelems);
DASSERT(x.nrows() == 1);

double DiagResult = 0, Result = 0;

// init xTemp and sum diag elements

for (int i = 0; i < nelems; i++)
	{
	xTemp[i] = x(i);
	DiagResult += A(i, i) * xTemp[i] * xTemp[i];
	}
// sum upper right triangle elems

for (i = 0; i < nelems; i++)
	{
	double xi = xTemp[i];
	for (int j = i+1; j < nelems; j++)
		Result += A(i, j) * xi * xTemp[j];
	}
Result *= 2;						// incorporate lower left triangle elements

return DiagResult + Result;
}

//-----------------------------------------------------------------------------
// This function is like DASSERT(A == A.t()) but tells you where the difference is, if any
// It's also quicker than testig A == A.t().

bool fIsMatrixSymmetric (const Mat &A, double Min)
{
bool fExact = (Min == 0.0);

if (A.nrows() != A.ncols())
	{
	lprintf("\nMatrix not square\n");
	return false;
	}
for (int iRow = 0; iRow < A.nrows(); iRow++)
	for (int iCol = iRow+1; iCol < A.ncols(); iCol++)
		if ((fExact && A(iRow, iCol) != A(iCol, iRow)) ||
		   (!fExact && !fEqual(A(iRow, iCol), A(iCol, iRow), Min)))
			{
			lprintf("\nMatrix not symmetric: A(%d,%d)=%.12g != A(%d,%d)=%.12g\n", iRow, iCol, A(iRow, iCol), iCol, iRow, A(iCol, iRow));
			//TODO PrintMat(A);
			return false;
			}

return true;
}

//-----------------------------------------------------------------------------
double SumElems (const Mat &m)
{
double sum = 0;
for (int iRow = 0; iRow < m.nrows(); iRow++)
	for (int iCol = 0; iCol < m.ncols(); iCol++)
		sum += m(iRow, iCol);
return sum;
}

//-----------------------------------------------------------------------------
double AbsSumElems (const Mat &m)
{
double sum = 0;
for (int iRow = 0; iRow < m.nrows(); iRow++)
	for (int iCol = 0; iCol < m.ncols(); iCol++)
		sum += fabs(m(iRow, iCol));
return sum;
}

}	// end namespace GslMat