Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Merge pull request #563 from karlnapf/master
some CMath convienience methods
  • Loading branch information
karlnapf committed Jun 3, 2012
2 parents 9b380e8 + 525eb66 commit f91331a
Show file tree
Hide file tree
Showing 5 changed files with 255 additions and 24 deletions.
99 changes: 76 additions & 23 deletions examples/undocumented/libshogun/mathematics_lapack.cpp
Expand Up @@ -23,8 +23,6 @@ bool is_equal(float64_t a, float64_t b, float64_t eps)

void test_ev()
{
#ifdef HAVE_LAPACK

SGMatrix<float64_t> A(3,3);
A(0,0)=0;
A(0,1)=1;
Expand Down Expand Up @@ -60,13 +58,69 @@ void test_ev()
ASSERT(is_equal(ev[0], -sqrt22*2, eps));
ASSERT(is_equal(ev[1], 0, eps));
ASSERT(is_equal(ev[2], sqrt22*2, eps));
}

#endif /* HAVE_LAPACK */
void test_matrix_multiply()
{
index_t n=10;
SGMatrix<float64_t> I=CMath::create_identity_matrix(n,1.0);

index_t m=4;
SGMatrix<float64_t> A(n, m);
CMath::range_fill_vector(A.matrix, m*n);
CMath::display_matrix(I, "I");
CMath::transpose_matrix(A.matrix, A.num_rows, A.num_cols);
CMath::display_matrix(A, "A transposed");
CMath::transpose_matrix(A.matrix, A.num_rows, A.num_cols);
CMath::display_matrix(A, "A");

SG_SPRINT("multiply A by I and check result\n");
SGMatrix<float64_t> A2=CMath::matrix_multiply(I, A);
ASSERT(A2.num_rows==A.num_rows);
ASSERT(A2.num_cols==A.num_cols);
CMath::display_matrix(A2);
for (index_t i=0; i<A2.num_rows; ++i)
{
for (index_t j=0; j<A2.num_cols; ++j)
ASSERT(A(i,j)==A2(i,j));
}

SG_SPRINT("multiply A by transposed I and check result\n");
SGMatrix<float64_t> A3=CMath::matrix_multiply(I, A, true);
ASSERT(A3.num_rows==I.num_rows);
ASSERT(A3.num_cols==A.num_cols);
CMath::display_matrix(A3);
for (index_t i=0; i<A2.num_rows; ++i)
{
for (index_t j=0; j<A2.num_cols; ++j)
ASSERT(A(i,j)==A3(i,j));
}

SG_SPRINT("multiply transposed A by I and check result\n");
SGMatrix<float64_t> A4=CMath::matrix_multiply(A, I, true, false);
ASSERT(A4.num_rows==A.num_cols);
ASSERT(A4.num_cols==I.num_cols);
CMath::display_matrix(A4);
for (index_t i=0; i<A.num_rows; ++i)
{
for (index_t j=0; j<A.num_cols; ++j)
ASSERT(A(i,j)==A4(j,i));
}

SG_SPRINT("multiply A by scaled I and check result\n");
SGMatrix<float64_t> A5=CMath::matrix_multiply(I, A, false, false, n);
ASSERT(A5.num_rows==I.num_rows);
ASSERT(A5.num_cols==A.num_cols);
CMath::display_matrix(A5);
for (index_t i=0; i<A2.num_rows; ++i)
{
for (index_t j=0; j<A2.num_cols; ++j)
ASSERT(n*A(i,j)==A5(i,j));
}
}

void test_lapack()
{
#ifdef HAVE_LAPACK
// size of square matrix
int N = 100;

Expand All @@ -93,10 +147,8 @@ void test_lapack()
status = 0;
wrap_dsygvx(1,'V','U',N,double_matrix,N,double_matrix,N,1,3,double_eigenvalues,double_eigenvectors,&status);
if (status!=0)
{
printf("DSYGVX/SSYGVX failed with code %d\n",status);
ASSERT(false);
}
SG_SERROR("DSYGVX/SSYGVX failed with code %d\n",status);

delete[] double_eigenvectors;

// DGEQRF+DORGQR
Expand All @@ -105,10 +157,8 @@ void test_lapack()
wrap_dgeqrf(N,N,double_matrix,N,double_tau,&status);
wrap_dorgqr(N,N,N,double_matrix,N,double_tau,&status);
if (status!=0)
{
printf("DGEQRF/DORGQR failed with code %d\n",status);
ASSERT(false);
}
SG_SERROR("DGEQRF/DORGQR failed with code %d\n",status);

delete[] double_tau;

// DGESVD
Expand All @@ -120,10 +170,8 @@ void test_lapack()
status = 0;
wrap_dgesvd('A','A',N,N,double_matrix,N,double_s,double_U,N,double_Vt,N,&status);
if (status!=0)
{
printf("DGESVD failed with code %d\n",status);
ASSERT(false);
}
SG_SERROR("DGESVD failed with code %d\n",status);

delete[] double_s;
delete[] double_U;
delete[] double_Vt;
Expand All @@ -132,24 +180,29 @@ void test_lapack()
status = 0;
wrap_dsyev('V','U',N,double_matrix,N,double_eigenvalues,&status);
if (status!=0)
{
printf("DSYEV failed with code %d\n",status);
ASSERT(false);
}
SG_SERROR("DSYEV failed with code %d\n",status);

delete[] double_eigenvalues;
delete[] double_matrix;

#endif // HAVE_LAPACK
}

int main(int argc, char** argv)
{
init_shogun();
init_shogun_with_defaults();

#ifdef HAVE_LAPACK
SG_SPRINT("checking lapack\n");
test_lapack();

SG_SPRINT("compute_eigenvectors\n");
test_ev();

SG_SPRINT("matrix_multiply\n");
test_matrix_multiply();
#endif

exit_shogun();
return 0;
}


120 changes: 119 additions & 1 deletion src/shogun/mathematics/Math.cpp
Expand Up @@ -164,6 +164,27 @@ void CMath::display_vector(const floatmax_t* vector, int32_t n,
SG_SPRINT("%s]\n", prefix);
}

template <>
void CMath::display_vector(const SGVector<float64_t> vector, const char* name,
const char* prefix)
{
CMath::display_vector(vector.vector, vector.vlen, name, prefix);
}

template <>
void CMath::display_vector(const SGVector<float32_t> vector, const char* name,
const char* prefix)
{
CMath::display_vector(vector.vector, vector.vlen, name, prefix);
}

template <>
void CMath::display_vector(const SGVector<int32_t> vector, const char* name,
const char* prefix)
{
CMath::display_vector(vector.vector, vector.vlen, name, prefix);
}

template <>
void CMath::display_matrix(
const int32_t* matrix, int32_t rows, int32_t cols, const char* name,
Expand Down Expand Up @@ -218,6 +239,72 @@ void CMath::display_matrix(
SG_SPRINT("%s]\n", prefix);
}

template <>
void CMath::display_matrix(
const SGMatrix<float64_t> matrix, const char* name,
const char* prefix)
{
CMath::display_matrix(matrix.matrix, matrix.num_rows, matrix.num_cols,
name, prefix);
}

template <>
void CMath::display_matrix(
const SGMatrix<float32_t> matrix, const char* name,
const char* prefix)
{
CMath::display_matrix(matrix.matrix, matrix.num_rows, matrix.num_cols,
name, prefix);
}

template <>
void CMath::display_matrix(
const SGMatrix<int32_t> matrix, const char* name,
const char* prefix)
{
CMath::display_matrix(matrix.matrix, matrix.num_rows, matrix.num_cols,
name, prefix);
}

template <>
SGMatrix<float64_t> CMath::create_identity_matrix(index_t size, float64_t scale)
{
SGMatrix<float64_t> I(size, size);
for (index_t i=0; i<size; ++i)
{
for (index_t j=0; j<size; ++j)
I(i,j)=i==j ? scale : 0.0;
}

return I;
}

template <>
SGMatrix<float32_t> CMath::create_identity_matrix(index_t size, float32_t scale)
{
SGMatrix<float32_t> I(size, size);
for (index_t i=0; i<size; ++i)
{
for (index_t j=0; j<size; ++j)
I(i,j)=i==j ? scale : 0.0;
}

return I;
}

template <>
SGMatrix<int32_t> CMath::create_identity_matrix(index_t size, int32_t scale)
{
SGMatrix<int32_t> I(size, size);
for (index_t i=0; i<size; ++i)
{
for (index_t j=0; j<size; ++j)
I(i,j)=i==j ? scale : 0;
}

return I;
}

}

SGVector<float64_t> CMath::fishers_exact_test_for_multiple_2x3_tables(SGMatrix<float64_t> tables)
Expand Down Expand Up @@ -488,7 +575,6 @@ float64_t CMath::entropy(float64_t* p, int32_t len)
return (float64_t) e;
}


//Howto construct the pseudo inverse (from "The Matrix Cookbook")
//
//Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
Expand Down Expand Up @@ -582,6 +668,38 @@ double* CMath::compute_eigenvectors(double* matrix, int n, int m)

return eigenvalues;
}

SGMatrix<float64_t> CMath::matrix_multiply(
SGMatrix<float64_t> A, SGMatrix<float64_t> B,
bool transpose_A, bool transpose_B, float64_t scale)
{
/* these variables store size of transposed matrices*/
index_t cols_A=transpose_A ? A.num_rows : A.num_cols;
index_t rows_A=transpose_A ? A.num_cols : A.num_rows;
index_t rows_B=transpose_B ? B.num_cols : B.num_rows;
index_t cols_B=transpose_B ? B.num_rows : B.num_cols;

/* do a dimension check */
if (cols_A!=rows_B)
{
SG_SERROR("CMath::matrix_multiply(): Dimension mismatch: "
"A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
}

/* allocate result matrix */
SGMatrix<float64_t> C(rows_A, cols_B);

/* multiply */
cblas_dgemm(CblasColMajor,
transpose_A ? CblasTrans : CblasNoTrans,
transpose_B ? CblasTrans : CblasNoTrans,
rows_A, cols_B, cols_A, scale,
A.matrix, A.num_rows, B.matrix, B.num_rows,
0.0, C.matrix, C.num_rows);

return C;
}

#endif //HAVE_LAPACK

float64_t CMath::dot(const float64_t* v1, const float64_t* v2, int32_t n)
Expand Down
29 changes: 29 additions & 0 deletions src/shogun/mathematics/Math.h
Expand Up @@ -493,6 +493,14 @@ class CMath : public CSGObject
}
}

/** returns the identity matrix, scaled by a factor
*
* @param size size of square idenity matrix
* @param scale (optional) scaling factor
*/
template <class T>
static SGMatrix<T> create_identity_matrix(index_t size, T scale);

#ifdef HAVE_LAPACK
/** compute eigenvalues and eigenvectors of symmetric matrix using
* LAPACK
Expand All @@ -514,6 +522,19 @@ class CMath : public CSGObject
* */
static double* compute_eigenvectors(double* matrix, int n, int m);

/* Computes scale* A*B, where A and B may be transposed.
* Asserts for matching inner dimensions.
* @param A matrix A
* @param transpose_A optional whether A should be transposed before
* @param B matrix B
* @param transpose_B optional whether B should be transposed before
* @param scale optional scaling factor for result
*/
static SGMatrix<float64_t> matrix_multiply(
SGMatrix<float64_t> A, SGMatrix<float64_t> B,
bool transpose_A=false, bool transpose_B=false,
float64_t scale=1.0);

/// inverses square matrix in-place
static void inverse(SGMatrix<float64_t> matrix);

Expand Down Expand Up @@ -1260,11 +1281,19 @@ class CMath : public CSGObject
const T* vector, int32_t n, const char* name="vector",
const char* prefix="");

template <class T> static void display_vector(
const SGVector<T>, const char* name="vector",
const char* prefix="");

/// display matrix (useful for debugging)
template <class T> static void display_matrix(
const T* matrix, int32_t rows, int32_t cols,
const char* name="matrix", const char* prefix="");

template <class T> static void display_matrix(
const SGMatrix<T> matrix, const char* name="matrix",
const char* prefix="");

/** performs a quicksort on an array output of length size
* it is sorted in ascending order
* (for type T1) and returns the index (type T2)
Expand Down
28 changes: 28 additions & 0 deletions src/shogun/statistics/QuadraticTimeMMD.cpp
Expand Up @@ -106,3 +106,31 @@ float64_t CQuadraticTimeMMD::compute_p_value(float64_t statistic)
return result;
}

SGVector<float64_t> CQuadraticTimeMMD::sample_null_spectrum(index_t num_samples)
{
SGVector<float64_t> null_samples(num_samples);

/* compute kernel matrix for XX, YY, XY */

/* imaginary matrix Kz=[K KL; KL' L] (MATLAB notation) */

/* centering matrix */

/* centering of matrix Kz=H*Kz*H */

/* compute eigenvalues */
/*kEigs = eigs(Kz,params.numEigs); %note: this retains only largest magnitude eigenvalues
%empirical eigenvalues scaled by 1/2/m: see p. 2 Shawe-Tayor
%et al. (2005)
kEigs = 1/2/m * abs(kEigs);
numEigs = length(kEigs); */

/* finally, sample from null distribution */
for (index_t i=0; i<num_samples; ++i)
{
/* 2*sum(kEigs.*(randn(length(kEigs),1)).^2); */
null_samples[i]=0;
}

return null_samples;
}

0 comments on commit f91331a

Please sign in to comment.