Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "DimensionalityReductionBase.h"
      23             : #include "core/ActionRegister.h"
      24             : 
      25             : //+PLUMEDOC DIMRED CLASSICAL_MDS
      26             : /*
      27             : Create a low-dimensional projection of a trajectory using the classical multidimensional
      28             :  scaling algorithm.
      29             : 
      30             : Multidimensional scaling (MDS) is similar to what is done when you make a map. You start with distances
      31             : between London, Belfast, Paris and Dublin and then you try to arrange points on a piece of paper so that the (suitably scaled)
      32             : distances between the points in your map representing each of those cities are related to the true distances between the cities.
      33             : Stating this more mathematically MDS endeavors to find an <a href="http://en.wikipedia.org/wiki/Isometry">isometry</a>
      34             : between points distributed in a high-dimensional space and a set of points distributed in a low-dimensional plane.
      35             : In other words, if we have \f$M\f$ \f$D\f$-dimensional points, \f$\mathbf{X}\f$,
      36             : and we can calculate dissimilarities between pairs them, \f$D_{ij}\f$, we can, with an MDS calculation, try to create \f$M\f$ projections,
      37             : \f$\mathbf{x}\f$, of the high dimensionality points in a \f$d\f$-dimensional linear space by trying to arrange the projections so that the
      38             : Euclidean distances between pairs of them, \f$d_{ij}\f$, resemble the dissimilarities between the high dimensional points.  In short we minimize:
      39             : 
      40             : \f[
      41             : \chi^2 = \sum_{i \ne j} \left( D_{ij} - d_{ij} \right)^2
      42             : \f]
      43             : 
      44             : where \f$D_{ij}\f$ is the distance between point \f$X^{i}\f$ and point \f$X^{j}\f$ and \f$d_{ij}\f$ is the distance between the projection
      45             : of \f$X^{i}\f$, \f$x^i\f$, and the projection of \f$X^{j}\f$, \f$x^j\f$.  A tutorial on this approach can be used to analyze simulations
      46             : can be found in the tutorial \ref lugano-5 and in the following <a href="https://www.youtube.com/watch?v=ofC2qz0_9_A&feature=youtu.be" > short video.</a>
      47             : 
      48             : \par Examples
      49             : 
      50             : The following command instructs plumed to construct a classical multidimensional scaling projection of a trajectory.
      51             : The RMSD distance between atoms 1-256 have moved is used to measure the distances in the high-dimensional space.
      52             : 
      53             : \plumedfile
      54             : data: COLLECT_FRAMES ATOMS=1-256
      55             : mat: EUCLIDEAN_DISSIMILARITIES USE_OUTPUT_DATA_FROM=data
      56             : mds: CLASSICAL_MDS USE_OUTPUT_DATA_FROM=mat NLOW_DIM=2
      57             : OUTPUT_ANALYSIS_DATA_TO_COLVAR USE_OUTPUT_DATA_FROM=mds FILE=rmsd-embed
      58             : \endplumedfile
      59             : 
      60             : The following section is for people who are interested in how this method works in detail. A solid understanding of this material is
      61             : not necessary to use MDS.
      62             : 
      63             : \section dim-sec Method of optimization
      64             : 
      65             : The stress function can be minimized using a standard optimization algorithm such as conjugate gradients or steepest descent.
      66             : However, it is more common to do this minimization using a technique known as classical scaling.  Classical scaling works by
      67             : recognizing that each of the distances $D_{ij}$ in the above sum can be written as:
      68             : 
      69             : \f[
      70             : D_{ij}^2 = \sum_{\alpha} (X^i_\alpha - X^j_\alpha)^2 = \sum_\alpha (X^i_\alpha)^2 + (X^j_\alpha)^2 - 2X^i_\alpha X^j_\alpha
      71             : \f]
      72             : 
      73             : We can use this expression and matrix algebra to calculate multiple distances at once.  For instance if we have three points,
      74             : \f$\mathbf{X}\f$, we can write distances between them as:
      75             : 
      76             : \f{eqnarray*}{
      77             : D^2(\mathbf{X}) &=& \left[ \begin{array}{ccc}
      78             : 0 & d_{12}^2 & d_{13}^2 \\
      79             : d_{12}^2 & 0 & d_{23}^2 \\
      80             : d_{13}^2 & d_{23}^2 & 0
      81             : \end{array}\right] \\
      82             : &=&
      83             : \sum_\alpha \left[ \begin{array}{ccc}
      84             : (X^1_\alpha)^2 & (X^1_\alpha)^2 & (X^1_\alpha)^2 \\
      85             : (X^2_\alpha)^2 & (X^2_\alpha)^2 & (X^2_\alpha)^2 \\
      86             : (X^3_\alpha)^2 & (X^3_\alpha)^2 & (X^3_\alpha)^2 \\
      87             : \end{array}\right]
      88             :  + \sum_\alpha \left[ \begin{array}{ccc}
      89             : (X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
      90             : (X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
      91             : (X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
      92             : \end{array}\right]
      93             : - 2 \sum_\alpha \left[ \begin{array}{ccc}
      94             : X^1_\alpha X^1_\alpha & X^1_\alpha X^2_\alpha & X^1_\alpha X^3_\alpha \\
      95             : X^2_\alpha X^1_\alpha & X^2_\alpha X^2_\alpha & X^2_\alpha X^3_\alpha \\
      96             : X^1_\alpha X^3_\alpha & X^3_\alpha X^2_\alpha & X^3_\alpha X^3_\alpha
      97             : \end{array}\right] \nonumber \\
      98             : &=& \mathbf{c 1^T} + \mathbf{1 c^T} - 2 \sum_\alpha \mathbf{x}_a \mathbf{x}^T_a =  \mathbf{c 1^T} + \mathbf{1 c^T} - 2\mathbf{X X^T}
      99             : \f}
     100             : 
     101             : This last equation can be extended to situations when we have more than three points.  In it \f$\mathbf{X}\f$ is a matrix that has
     102             : one high-dimensional point on each of its rows and \f$\mathbf{X^T}\f$ is its transpose.  \f$\mathbf{1}\f$ is an \f$M \times 1\f$ vector
     103             : of ones and \f$\mathbf{c}\f$ is a vector with components given by:
     104             : 
     105             : \f[
     106             : c_i = \sum_\alpha (x_\alpha^i)^2
     107             : \f]
     108             : 
     109             : These quantities are the diagonal elements of \f$\mathbf{X X^T}\f$, which is a dot product or Gram Matrix that contains the
     110             : dot product of the vector \f$X_i\f$ with the vector \f$X_j\f$ in element \f$i,j\f$.
     111             : 
     112             : In classical scaling we introduce a centering matrix \f$\mathbf{J}\f$ that is given by:
     113             : 
     114             : \f[
     115             : \mathbf{J} = \mathbf{I} - \frac{1}{M} \mathbf{11^T}
     116             : \f]
     117             : 
     118             : where \f$\mathbf{I}\f$ is the identity.  Multiplying the equations above from the front and back by this matrix and a factor of a \f$-\frac{1}{2}\f$ gives:
     119             : 
     120             : \f{eqnarray*}{
     121             :  -\frac{1}{2} \mathbf{J} \mathbf{D}^2(\mathbf{X}) \mathbf{J} &=& -\frac{1}{2}\mathbf{J}( \mathbf{c 1^T} + \mathbf{1 c^T} - 2\mathbf{X X^T})\mathbf{J} \\
     122             :  &=& -\frac{1}{2}\mathbf{J c 1^T J} - \frac{1}{2} \mathbf{J 1 c^T J} + \frac{1}{2} \mathbf{J}(2\mathbf{X X^T})\mathbf{J} \\
     123             :  &=& \mathbf{ J X X^T J } = \mathbf{X X^T } \label{eqn:scaling}
     124             : \f}
     125             : 
     126             : The fist two terms in this expression disappear because \f$\mathbf{1^T J}=\mathbf{J 1} =\mathbf{0}\f$, where \f$\mathbf{0}\f$
     127             : is a matrix containing all zeros.  In the final step meanwhile we use the fact that the matrix of squared distances will not
     128             : change when we translate all the points.  We can thus assume that the mean value, \f$\mu\f$, for each of the components, \f$\alpha\f$:
     129             : \f[
     130             : \mu_\alpha = \frac{1}{M} \sum_{i=1}^N \mathbf{X}^i_\alpha
     131             : \f]
     132             : is equal to 0 so the columns of \f$\mathbf{X}\f$ add up to 0.  This in turn means that each of the columns of
     133             : \f$\mathbf{X X^T}\f$ adds up to zero, which is what allows us to write \f$\mathbf{ J X X^T J } = \mathbf{X X^T }\f$.
     134             : 
     135             : The matrix of squared distances is symmetric and positive-definite we can thus use the spectral decomposition to decompose it as:
     136             : 
     137             : \f[
     138             : \Phi= \mathbf{V} \Lambda \mathbf{V}^T
     139             : \f]
     140             : 
     141             : Furthermore, because the matrix we are diagonalizing, \f$\mathbf{X X^T}\f$, is the product of a matrix and its transpose
     142             : we can use this decomposition to write:
     143             : 
     144             : \f[
     145             : \mathbf{X} =\mathbf{V} \Lambda^\frac{1}{2}
     146             : \f]
     147             : 
     148             : Much as in PCA there are generally a small number of large eigenvalues in \f$\Lambda\f$ and many small eigenvalues.
     149             : We can safely use only the large eigenvalues and their corresponding eigenvectors to express the relationship between
     150             : the coordinates \f$\mathbf{X}\f$.  This gives us our set of low-dimensional projections.
     151             : 
     152             : This derivation makes a number of assumptions about the how the low dimensional points should best be arranged to minimize
     153             : the stress. If you use an interactive optimization algorithm such as SMACOF you may thus be able to find a better
     154             : (lower-stress) projection of the points.  For more details on the assumptions made
     155             : see <a href="http://quest4rigor.com/tag/multidimensional-scaling/"> this website.</a>
     156             : */
     157             : //+ENDPLUMEDOC
     158             : 
     159             : namespace PLMD {
     160             : namespace dimred {
     161             : 
     162             : class ClassicalMultiDimensionalScaling : public DimensionalityReductionBase {
     163             : public:
     164             :   static void registerKeywords( Keywords& keys );
     165             :   explicit ClassicalMultiDimensionalScaling( const ActionOptions& ao );
     166             :   void calculateProjections( const Matrix<double>&, Matrix<double>& ) override;
     167             : };
     168             : 
     169       10433 : PLUMED_REGISTER_ACTION(ClassicalMultiDimensionalScaling,"CLASSICAL_MDS")
     170             : 
     171           8 : void ClassicalMultiDimensionalScaling::registerKeywords( Keywords& keys ) {
     172           8 :   DimensionalityReductionBase::registerKeywords( keys );
     173           8 : }
     174             : 
     175           7 : ClassicalMultiDimensionalScaling::ClassicalMultiDimensionalScaling( const ActionOptions& ao):
     176             :   Action(ao),
     177           7 :   DimensionalityReductionBase(ao)
     178             : {
     179           7 :   if( dimredbase ) error("input to CLASSICAL_MDS should not be output from dimensionality reduction object");
     180           7 : }
     181             : 
     182          12 : void ClassicalMultiDimensionalScaling::calculateProjections( const Matrix<double>& targets, Matrix<double>& projections ) {
     183             :   // Retrieve the distances from the dimensionality reduction object
     184          12 :   double half=(-0.5); Matrix<double> distances( half*targets );
     185             : 
     186             :   // Apply centering transtion
     187             :   unsigned n=distances.nrows(); double sum;
     188             :   // First HM
     189        1564 :   for(unsigned i=0; i<n; ++i) {
     190      389056 :     sum=0; for(unsigned j=0; j<n; ++j) sum+=distances(i,j);
     191      389056 :     for(unsigned j=0; j<n; ++j) distances(i,j) -= sum/n;
     192             :   }
     193             :   // Now (HM)H
     194        1564 :   for(unsigned i=0; i<n; ++i) {
     195      389056 :     sum=0; for(unsigned j=0; j<n; ++j) sum+=distances(j,i);
     196      389056 :     for(unsigned j=0; j<n; ++j) distances(j,i) -= sum/n;
     197             :   }
     198             : 
     199             :   // Diagonalize matrix
     200          12 :   std::vector<double> eigval(n); Matrix<double> eigvec(n,n);
     201          12 :   diagMat( distances, eigval, eigvec );
     202             : 
     203             :   // Pass final projections to map object
     204        1564 :   for(unsigned i=0; i<n; ++i) {
     205        4654 :     for(unsigned j=0; j<projections.ncols(); ++j) projections(i,j)=std::sqrt(eigval[n-1-j])*eigvec(n-1-j,i);
     206             :   }
     207          12 : }
     208             : 
     209             : }
     210             : }