72 assert( dimensions >= 1 );
73 assert( iterations >= 1 );
74 assert( initial_values.
rows() == distances.
rows() );
75 assert( initial_values.
cols() == dimensions );
78 double const learning_rate = 0.05;
79 double const r_metric = 2.0;
80 size_t const n = distances.
rows();
83 auto result = initial_values;
84 assert( result.rows() == n );
85 assert( result.cols() == dimensions );
90 auto dhdum = std::vector<double>( n, 0.0 );
95 for(
size_t j = 0; j < perms.
cols(); ++j ) {
96 std::iota( perms.
col(j).begin(), perms.
col(j).end(), 0 );
101 for(
size_t it = 0; it < iterations; ++it ) {
104 for(
size_t rp = 0; rp < n; ++rp ) {
105 auto const m = perms( rp, it );
108 for(
size_t i = 0; i < n; ++i ) {
109 for(
size_t j = 0; j < dimensions; ++j ) {
110 pmat( i, j ) = result( m, j ) - result( i, j );
113 for(
size_t i = 0; i < n; ++i ) {
115 for(
size_t j = 0; j < dimensions; ++j ) {
116 temp += std::pow( std::abs( pmat( i, j )), r_metric );
118 dhdum[ i ] = std::pow( temp, 1.0 / r_metric );
120 for(
size_t i = 0; i < n; ++i ) {
124 dh( m, i ) = dhdum[ i ];
125 dh( i, m ) = dhdum[ i ];
127 for(
size_t i = 0; i < n - 1; ++i ) {
128 size_t const ii = ( i < m ) ? i : i + 1;
131 * ( dhdum[ ii ] - distances( ii, m ) )
132 * std::pow( dhdum[ ii ], 1.0 - r_metric )
134 for(
size_t j = 0; j < dimensions; ++j ) {
135 dhmat( i, j ) = temp;
138 for(
size_t i = 0; i < n - 1; ++i ) {
139 size_t const ii = ( i < m ) ? i : i + 1;
140 for(
size_t j = 0; j < dimensions; ++j ) {
141 double temp = result( ii, j );
142 temp += dhmat( i, j )
143 * std::pow( std::abs( pmat( ii, j )), r_metric - 1.0 )
144 * signum<double>( pmat( ii, j ))
146 result( ii, j ) = temp;
164 assert( dimensions >= 1 );
165 assert( iterations >= 1 );
166 assert( initial_values.
rows() == distances.
rows() );
167 assert( initial_values.
cols() == dimensions );
171 auto result = initial_values;
172 auto previous = initial_values;
176 assert( idists.rows() == distances.
rows() );
177 assert( idists.cols() == distances.
cols() );
183 for(
size_t it = 0; it < iterations; ++it ) {
186 for(
size_t i = 0; i < distances.
rows(); ++i ) {
187 for(
size_t j = 0; j < distances.
cols(); ++j ) {
188 if( i == j || std::abs( idists(i,j)) <
MDS_EPSILON ) {
189 stress( i, j ) = 0.0;
191 stress( i, j ) = - distances( i, j ) / idists( i, j );
197 for(
size_t j = 0; j < distances.
cols(); ++j ) {
199 for(
size_t i = 0; i < distances.
rows(); ++i ) {
200 temp += stress( i, j );
202 stress( j, j ) = - temp;
206 assert( stress.rows() == result.rows() );
207 assert( previous.rows() == stress.cols() );
208 assert( previous.cols() == result.cols() );
209 for(
size_t i = 0; i < result.rows(); ++i ) {
210 for(
size_t j = 0; j < result.cols(); ++j ) {
212 for(
size_t k = 0; k < stress.cols(); ++k ) {
213 temp += stress( i, k ) * previous( k, j );
215 result( i, j ) = temp /
static_cast<double>( distances.
rows() );
243 auto distrib = std::uniform_real_distribution<double>( 0.0, 1.0 );
245 for(
auto& e : initial ) {
246 auto const r = distrib( engine );
250 mean /= initial.size();
251 assert( 0.0 <= mean && mean <= 1.0 );
254 for(
auto& e : initial ) {
255 e *= 0.1 * mean / ( 1.0 / 3.0 * std::sqrt( static_cast<double>( dimensions )));
270 throw std::invalid_argument(
"MDS input distance matrix is not square." );
272 if( dimensions < 1 ) {
273 throw std::invalid_argument(
"MDS dimensions has to be >= 1." );
275 if( iterations < 1 ) {
276 throw std::invalid_argument(
"MDS number of iterations has to be >= 1." );
278 if( initial_values.
rows() != distances.
rows() || initial_values.
cols() != dimensions ) {
279 throw std::invalid_argument(
"MDS initial values matrix has invalid dimensions." );
281 if( distances.
empty() ) {
285 switch( algorithm ) {
288 distances, initial_values, dimensions, iterations
293 distances, initial_values, dimensions, iterations
Provides some valuable algorithms that are not part of the C++ 11 STL.
static Matrix< double > multi_dimensional_scaling_ucf(Matrix< double > const &distances, Matrix< double > const &initial_values, size_t dimensions, size_t iterations)
Use the UCF implementation (recommended).
constexpr double MDS_EPSILON
Container namespace for all symbols of genesis in order to keep them separate when used as a library...
bool is_square(Matrix< T > const &data)
Return whether a Matrix is a square matrix, that is, whether its number of rows and number of columns...
static Matrix< double > multi_dimensional_scaling_smacof(Matrix< double > const &distances, Matrix< double > const &initial_values, size_t dimensions, size_t iterations)
double mean(const Histogram &h)
Compute the bin-weighted arithmetic mean.
MatrixCol< self_type, value_type > col(size_t col)
MdsAlgorithm
Choice of algorithm to use for Multi-Dimensional Scaling (MDS).
Matrix< double > euclidean_distance_matrix(Matrix< double > const &data)
Calculate the pairwise euclidean distance matrix between the rows of a given matrix.
Matrix< double > multi_dimensional_scaling(Matrix< double > const &distances, size_t dimensions, size_t iterations, MdsAlgorithm algorithm)
Multi-Dimensional Scaling (MDS).
std::default_random_engine & random_engine()
Returns the default engine for random number generation.
Use the SMACOF implementation.
static Options & get()
Returns a single instance of this class.