Permalink
Browse files

added sparse-matrix linear operator in log-det framework

  • Loading branch information...
lambday committed Jul 10, 2013
1 parent 72301ea commit d393aaf5b5004ee6926fdf49d74457ef9a1bbe43
@@ -3986,7 +3986,7 @@ bool TParameter::copy(TParameter* target)
else
num_bytes=*m_datatype.m_length_y *
(*m_datatype.m_length_y) * m_datatype.sizeof_stype();
*(char**)target->m_parameter=SG_MALLOC(char, num_bytes);
*(char**)target->m_parameter=SG_CALLOC(char, num_bytes);
/* use length of source */
*target->m_datatype.m_length_y=*m_datatype.m_length_y;
@@ -66,6 +66,20 @@ template <class T> class SGSparseMatrix : public SGReferencedData
return *this;
}
/* compute sparse-matrix dense-vector multiplication
* @param v the dense-vector to be multiplied with
* @return the result vector \f$Q*v\f$, Q being this sparse matrix
*/
const SGVector<T> operator*(SGVector<T> v) const
{
SGVector<T> result(num_vectors);
ASSERT(v.vlen==num_features);
for (index_t i=0; i<num_vectors; ++i)
result[i]=sparse_matrix[i].dense_dot(1.0, v.vector, v.vlen, 0.0);
return result;
}
/** load sparse matrix from file
*
* @param loader File object via which to load data
@@ -113,7 +113,7 @@ SGVector<T> CDenseMatrixOperator<T>::apply(SGVector<T> b) const
Map<MatrixXt> _op(m_operator.matrix, m_operator.num_rows,
m_operator.num_cols);
SGVector<T> result(b.vlen);
SGVector<T> result(m_operator.num_rows);
Map<VectorXt> _result(result.vector, result.vlen);
_result=_op*_b;
@@ -0,0 +1,171 @@
/*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* Written (W) 2013 Soumyajit De
*/
#include <shogun/lib/config.h>
#include <shogun/lib/SGVector.h>
#include <shogun/lib/SGSparseMatrix.h>
#include <shogun/lib/SGSparseVector.h>
#include <shogun/base/Parameter.h>
#include <shogun/mathematics/logdet/SparseMatrixOperator.h>
namespace shogun
{
template<class T>
CSparseMatrixOperator<T>::CSparseMatrixOperator()
: CMatrixOperator<T>()
{
init();
SG_SGCDEBUG("%s created (%p)\n", this->get_name(), this);
}
template<class T>
CSparseMatrixOperator<T>::CSparseMatrixOperator(SGSparseMatrix<T> op)
: CMatrixOperator<T>(op.num_features),
m_operator(op)
{
init();
SG_SGCDEBUG("%s created (%p)\n", this->get_name(), this);
}
template<class T>
void CSparseMatrixOperator<T>::init()
{
CSGObject::set_generic<T>();
CSGObject::m_parameters->add(&m_operator,
"sparse_matrix",
"The sparse matrix of the linear operator");
}
template<class T>
CSparseMatrixOperator<T>::~CSparseMatrixOperator()
{
SG_SGCDEBUG("%s destroyed (%p)\n", this->get_name(), this);
}
template<class T>
SGSparseMatrix<T> CSparseMatrixOperator<T>::get_matrix_operator() const
{
return m_operator;
}
template<class T>
SGVector<T> CSparseMatrixOperator<T>::get_diagonal() const
{
REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
const int32_t diag_size=m_operator.num_vectors>m_operator.num_features ?
m_operator.num_features : m_operator.num_vectors;
SGVector<T> diag(diag_size);
diag.set_const(static_cast<T>(0));
for (index_t i=0; i<diag_size; ++i)
{
SGSparseVectorEntry<T>* current_row=m_operator[i].features;
for (index_t j=0; j<m_operator[i].num_feat_entries; ++j)
{
if (i==current_row[j].feat_index)
{
diag[i]=current_row[j].entry;
break;
}
}
}
return diag;
}
template<class T>
void CSparseMatrixOperator<T>::set_diagonal(SGVector<T> diag)
{
REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
REQUIRE(diag.vector, "Diagonal not initialized!\n");
const int32_t diag_size=m_operator.num_vectors>m_operator.num_features ?
m_operator.num_features : m_operator.num_vectors;
REQUIRE(diag_size==diag.vlen, "Dimension mismatch!\n");
bool need_sorting=false;
for (index_t i=0; i<diag_size; ++i)
{
SGSparseVectorEntry<T>* current_row=m_operator[i].features;
bool inserted=false;
// we just change the entry if the diagonal element for this row exists
for (index_t j=0; j<m_operator[i].num_feat_entries; ++j)
{
if (i==current_row[j].feat_index)
{
current_row[j].entry=diag[i];
inserted=true;
break;
}
}
// we create a new entry if the diagonal element for this row doesn't exist
if (!inserted)
{
index_t j=m_operator[i].num_feat_entries;
m_operator[i].num_feat_entries=j+1;
m_operator[i].features=SG_REALLOC(SGSparseVectorEntry<T>,
m_operator[i].features, j, j+1);
m_operator[i].features[j].feat_index=i;
m_operator[i].features[j].entry=diag[i];
need_sorting=true;
}
}
if (need_sorting)
m_operator.sort_features();
}
template<class T>
SGVector<T> CSparseMatrixOperator<T>::apply(SGVector<T> b) const
{
REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
REQUIRE(CLinearOperator<T>::m_dimension==b.vlen,
"Number of rows of vector must be equal to the "
"number of cols of the operator!\n");
SGVector<T> result(m_operator.num_vectors);
result=m_operator*b;
return result;
}
#define UNDEFINED(type) \
template<> \
SGVector<type> CSparseMatrixOperator<type>::apply(SGVector<type> b) const \
{ \
SG_SERROR("Not supported for %s\n", #type);\
return b; \
}
UNDEFINED(bool)
UNDEFINED(char)
#undef UNDEFINED
template class CSparseMatrixOperator<bool>;
template class CSparseMatrixOperator<char>;
template class CSparseMatrixOperator<int8_t>;
template class CSparseMatrixOperator<uint8_t>;
template class CSparseMatrixOperator<int16_t>;
template class CSparseMatrixOperator<uint16_t>;
template class CSparseMatrixOperator<int32_t>;
template class CSparseMatrixOperator<uint32_t>;
template class CSparseMatrixOperator<int64_t>;
template class CSparseMatrixOperator<uint64_t>;
template class CSparseMatrixOperator<float32_t>;
template class CSparseMatrixOperator<float64_t>;
template class CSparseMatrixOperator<floatmax_t>;
template class CSparseMatrixOperator<complex64_t>;
}
@@ -0,0 +1,87 @@
/*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* Written (W) 2013 Soumyajit De
*/
#ifndef SPARSE_MATRIX_OPERATOR_H_
#define SPARSE_MATRIX_OPERATOR_H_
#include <shogun/lib/config.h>
#include <shogun/mathematics/logdet/MatrixOperator.h>
namespace shogun
{
template<class T> class SGVector;
template<class T> class SGSparseMatrix;
/** @brief Class that represents a sparse-matrix linear operator.
* It computes matrix-vector product \f$Ax\f$ in its apply method,
* \f$A\in\mathbb{C}^{m\times n},A:\mathbb{C}^{n}\rightarrow \mathbb{C}^{m}\f$
* being the matrix operator and \f$x\in\mathbb{C}^{n}\f$ being the vector.
* The result is a vector \f$y\in\mathbb{C}^{m}\f$.
*/
template<class T> class CSparseMatrixOperator : public CMatrixOperator<T>
{
/** this class has support for complex64_t */
typedef bool supports_complex64_t;
public:
/** default constructor */
CSparseMatrixOperator();
/**
* constructor
*
* @param op the sparse matrix to be used as the linear operator
*/
CSparseMatrixOperator(SGSparseMatrix<T> op);
/** destructor */
~CSparseMatrixOperator();
/**
* method that applies the sparse-matrix linear operator to a vector
*
* @param b the vector to which the linear operator applies
* @return the result vector
*/
virtual SGVector<T> apply(SGVector<T> b) const;
/**
* method that sets the main diagonal of the matrix
*
* @param diag the diagonal to be set
*/
virtual void set_diagonal(SGVector<T> diag);
/**
* method that returns the main diagonal of the matrix
*
* @return the diagonal
*/
virtual SGVector<T> get_diagonal() const;
/** @return the sparse matrix operator */
SGSparseMatrix<T> get_matrix_operator() const;
/** @return object name */
virtual const char* get_name() const
{
return "SparseMatrixOperator";
}
private:
/** the sparse matrix operator */
SGSparseMatrix<T> m_operator;
/** initialize with default values and register params */
void init();
};
}
#endif // SPARSE_MATRIX_OPERATOR_H_
Oops, something went wrong.

0 comments on commit d393aaf

Please sign in to comment.