GenMat.h

Go to the documentation of this file.

// -*- C++ -*-

// PLearn (A C++ Machine Learning Library)
// Copyright (C) 1998 Pascal Vincent
// Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal
//

// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
// 
//  1. Redistributions of source code must retain the above copyright
//     notice, this list of conditions and the following disclaimer.
// 
//  2. Redistributions in binary form must reproduce the above copyright
//     notice, this list of conditions and the following disclaimer in the
//     documentation and/or other materials provided with the distribution.
// 
//  3. The name of the authors may not be used to endorse or promote
//     products derived from this software without specific prior written
//     permission.
// 
// THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
// NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
// 
// This file is part of the PLearn library. For more information on the PLearn
// library, go to the PLearn Web site at www.plearn.org


 

/* *******************************************************      
   * $Id: GenMat.h,v 1.2 2004/02/20 21:11:46 chrish42 Exp $
   * Generic Matrix (template) operations
   * AUTHORS: Yoshua Bengio
   * This file is part of the PLearn library.
   ******************************************************* */


#ifndef GenMat_INC
#define GenMat_INC

namespace PLearn {
using namespace std;


template<class MatT>
class SquaredSymmMatT 
{
protected:
  MatT& A_;
  Vec Ax;
public:
  SquaredSymmMatT(MatT& A) : A_(A), Ax(A.length()) 
  { 
    if (A.length() != A.width()) 
      PLERROR("SquaredSymmMatT expects a square symmetric matrix"); 
  }
  int length() const { return A_.length(); }
  int width() const { return A_.width(); }
  void product(const Vec& x, Vec& y)
  {
    product(A_, x,Ax);
    product(A_, Ax,y);
  }

  void diag(Vec& d) { diagonalOfSquare(A_, d); }
  void diagonalOfSquare(Vec& d) 
  { PLERROR("SquaredSymmMatT::diagonalOfSquare not implemented"); }

};

template<class MatT>
class ReverseMatT 
{
protected:
  MatT& A_;
  real alpha_;
public:
  ReverseMatT(MatT& A, real alpha) : A_(A), alpha_(alpha)
  { 
    if (A.length() != A.width()) 
      PLERROR("ReverseMatT expects a square symmetric matrix"); 
  }
  int length() const { return A_.length(); }
  int width() const { return A_.width(); }
  void product(const Vec& x, Vec& y)
  {
    product(A_, x,y);
    y*=-1;
    multiplyAcc(y, x,alpha_);
  }

  void diag(Vec& d) 
  { 
    diag(A_, d); 
    diag*=-1;
    diag+=alpha_;
  }
  void diagonalOfSquare(Vec& d) 
  {
    Vec diag_A(A_.length());
    diag(A_, diag_A);
    diagonalOfSquare(A_, d);
    multiplyAcc(d, diag_A,-2*alpha_);
    d+=alpha_*alpha_;
  }

};

template<class MatT>
class MatTPlusSumSquaredVec
{
public:
  MatT& A_;
  Mat X_; 
public:
  MatTPlusSumSquaredVec(MatT& A, int alloc_n_vectors) : A_(A), X_(alloc_n_vectors,A.length())
  { X_.resize(0,A.length()); }
  MatTPlusSumSquaredVec(MatT& A, Mat& X) : A_(A), X_(X) 
  { 
    if (X.width()!=A.length() || A.length()!=A.width()) 
      PLERROR("MatTPlusSumSquaredVec(MatT,Mat) arguments don't match"); 
  }
  void squaredVecAcc(Vec& x)
  {
    X_.resize(X_.length()+1,X_.width());
    X_(X_.length()-1) << x;
  }
  int length() const { return A_.length(); }
  int width() const { return A_.width(); }
  void product(const Vec& x, Vec& y)
  {
    product(A_, x,y);
    for (int t=0;t<X_.length();t++)
    {
      Vec x_t = X_(t);
      multiplyAcc(y, x_t,dot(x_t,x));
    }
  }

  void diag(Vec& d) 
  { 
    diag(A_, d); 
    for (int t=0;t<X_.length();t++)
      squareAcc(d, X_(t));
  }
  void diagonalOfSquare(Vec& d) 
  { PLERROR("MatTPlusSumSquaredVec::diagonalOfSquare not implemented"); }
};


#if 0
// for debugging
#define MatT Mat
inline bool SolveLinearSymmSystemByCG(MatT A, Vec x, const Vec& b, int& max_iter, real& tol, real lambda)
{
  real resid;
  int n=A.nrows();
  Vec p(n), z(n), q(n), invdiag(n);
  real alpha, beta, rho, previous_rho;
  real normb = norm(b);

  // inverse diagonal of A, for preconditionning
  diag(A, invdiag);
  if (lambda>0)
    invdiag+=lambda;
  invertElements(invdiag);

  Vec r(n);
  // r = b - (A+lambda*I)*x;
  product(A, x,r);
  if (lambda>0)
    multiplyAcc(r, x,lambda);
  r*=-1;
  r+=b;

  if (normb == 0.0) 
    normb = 1;
  
  resid = norm(r);
  //cout << "at 0, |r| = " << resid << endl;
  if ((resid / normb) <= tol) {
    tol = resid;
    max_iter = 0;
    return true;
  }

  for (int i = 1; i <= max_iter; i++) {
    // z = M.solve(r); i.e. solve M z = r, i.e. diag(A+lambda I) z = r
    // i.e. z_i = r_i / (A_{i,i} + lambda)
    multiply(r,invdiag,z);

    rho = dot(r, z);
    
    if (i == 1)
      p << z;
    else {
      beta = rho / previous_rho;
      // p = z + beta * p;
      multiplyAdd(z,p,beta,p);
    }
    
    // q = (A+lambda I)*p;
    product(A, p,q);
    multiplyAcc(q, p,lambda);

    alpha = rho / dot(p, q);

    // x += alpha * p;
    multiplyAcc(x, p,alpha);
    // r -= alpha * q;
    multiplyAcc(r, q,-alpha);
    resid = norm(r);
    // cout << "at " << i << ", |r| = " << resid << endl;
    if ((resid / normb) <= tol) {
      tol = resid;
      max_iter = i;
      return true;     
    }

    previous_rho = rho;
  }
  
  tol = resid;
  return false;
}
#undef MatT
#else

template <class MatT>
bool SolveLinearSymmSystemByCG(MatT A, Vec x, const Vec& b, int& max_iter, real& tol, real lambda)
{
  real resid;
  int n=A.length();
  Vec p(n), z(n), q(n), invdiag(n);
  real alpha, beta, rho, previous_rho;
  real normb = norm(b);

  diag(A, invdiag);
  if (lambda>0)
    invdiag+=lambda;
  invertElements(invdiag);

  Vec r(n);
  product(A, x,r);
  if (lambda>0)
    multiplyAcc(r, x,lambda);
  r*=-1;
  r+=b;

  if (normb == 0.0) 
    normb = 1;
  
  resid = norm(r);
  //cout << "at 0, |r| = " << resid << endl;
  if ((resid / normb) <= tol) {
    tol = resid;
    max_iter = 0;
    return true;
  }

  for (int i = 1; i <= max_iter; i++) {
    multiply(r,invdiag,z);

    rho = dot(r, z);
    
    if (i == 1)
      p << z;
    else {
      beta = rho / previous_rho;
      multiplyAdd(z,p,beta,p);
    }
    
    product(A, p,q);
    multiplyAcc(q, p,lambda);

    alpha = rho / dot(p, q);

    multiplyAcc(x, p,alpha);
    multiplyAcc(r, q,-alpha);
    resid = norm(r);
    if ((resid / normb) <= tol) {
      tol = resid;
      max_iter = i;
      return true;     
    }

    previous_rho = rho;
  }
  
  tol = resid;
  return false;
}
#endif

#if 0

#define MatT Mat
#if 0
#define MatT Mat
inline real PowerIteration(MatT& A, Vec x0, int& n_iterations, 
                    real RayleighQuotientTolerance, Mat final_vectors, 
                    int& final_offset)
{
  int N=final_vectors.length();
  if (N<3) PLERROR("PowerIteration: final_vectors.length_ = %d should be >= 3",N);
  Vec previous=final_vectors(0);
  Vec current=final_vectors(1);
  Vec next=final_vectors(2);
  previous << x0; 
  product(A, previous,current);
  real current_norm=norm(current), next_norm, current_Rq, previous_Rq=0;
  current/=current_norm;
  for (int it=1;it<=n_iterations;it++)
  {
    product(A, current,next);
    // check Rayleigh quotient (note that norm(current)=1)
    current_Rq = dot(current,next);
    // normalize
    next_norm = norm(next);
    next/=next_norm;
    cout << "at iteration " << it << ", R(A,x) = " << current_Rq << ", |Ax|/|x| = " 
         << next_norm << endl;
    if (current_Rq < previous_Rq)
      PLWARNING("PowerIteration: something strange, x'Ax/x'x is decreasing %g->%g",
              previous_Rq, current_Rq);
    if (it>=N && current_Rq / previous_Rq - 1 < RayleighQuotientTolerance)
    {
      // stop
      n_iterations = it;
      final_offset = it%N;
      return current_Rq;
    }
    // iterate
    previous_Rq = current_Rq;
    current_norm = next_norm;
    previous = current;
    current = next;
    next = final_vectors((it+2)%N);
  }
}
#undef MatT
#else

template<class MatT>
real PowerIteration(MatT& A, Vec x0, int& n_iterations, 
                    real RayleighQuotientTolerance, Mat final_vectors, 
                    int& final_offset)
{
  int N=final_vectors.length();
  if (N<3) PLERROR("PowerIteration: final_vectors.length_ = %d should be >= 3",N);
  Vec previous=final_vectors(0);
  Vec current=final_vectors(1);
  Vec next=final_vectors(2);
  previous << x0; 
  product(A, previous,current);
  real current_norm=norm(current), next_norm, current_Rq, previous_Rq=0;
  current/=current_norm;
  for (int it=1;it<=n_iterations;it++)
  {
    product(A, current,next);
    current_Rq = dot(current,next);
    next_norm = norm(next);
    next/=next_norm;
    //cout << "at iteration " << it << ", R(A,x) = " << current_Rq << ", |Ax|/|x| = " 
    if (current_Rq < previous_Rq)
      PLWARNING("PowerIteration: something strange, x'Ax/x'x is decreasing %g->%g",
              previous_Rq, current_Rq);
    if (it>=N && current_Rq / previous_Rq - 1 < RayleighQuotientTolerance)
    {
      n_iterations = it;
      final_offset = it%N;
      return next_norm;
    }
    previous_Rq = current_Rq;
    current_norm = next_norm;
    previous = current;
    current = next;
    next = final_vectors((it+2)%N);
  }
}
#endif

  int GramSchmidtOrthogonalization(Mat A, real tolerance=1e-6);

template <class MatT>
real PowerIteration(MatT A, Vec x0, int& n_iterations, 
                    real RayleighQuotientTolerance, Mat final_vectors, 
                    int& final_offset, bool verbose=false)
{
  int N=final_vectors.length();
  if (N<3) PLERROR("PowerIteration: final_vectors.length_ = %d should be >= 3",N);
  Vec previous=final_vectors(0);
  Vec current=final_vectors(1);
  Vec next=final_vectors(2);
  previous << x0; 
  real max_x = max(previous);
  if (max_x<0) previous*=-1;
  product(A, previous,current);
  real current_norm=norm(current), next_norm, current_Rq, previous_Rq=0;
  max_x = max(current);
  if (max_x<0) current*=-1;
  current/=current_norm;
  int it=1;
  for (;it<=n_iterations;it++)
  {
    product(A, current,next);
    current_Rq = dot(current,next);
    next_norm = norm(next);
    next/=next_norm;
    max_x = max(next);
    if (max_x<0) next*=-1;
    if (verbose)
    {
      cout << "at iteration " << it << ", R(A,x) = " << current_Rq << ", |Ax|/|x| = " 
           << next_norm << endl;
    }
    if (current_Rq < previous_Rq)
      PLWARNING("PowerIteration: something strange, x'Ax/x'x is decreasing %g->%g",
              previous_Rq, current_Rq);
    if (it>=N && current_Rq / previous_Rq - 1 < RayleighQuotientTolerance)
    {
      n_iterations = it;
      final_offset = it%N;
      if (verbose)
        cout << "power iteration finishes with |Ax|/|x| = " << next_norm << endl;
      return next_norm;
    }
    previous_Rq = current_Rq;
    current_norm = next_norm;
    previous = current;
    current = next;
    next = final_vectors((it+2)%N);
  }
  final_offset = it%N;
  if (verbose)
    cout << "power iteration finishes FAILING with |Ax|/|x| = " << next_norm << endl;
  return next_norm;
}

template<class MatT>
real InversePowerIteration(MatT A, Vec x0, int& n_iterations, 
                           int max_n_cg_iter,
                           real RTolerance, Mat final_vectors, 
                           int& final_offset, real regularizer, bool verbose=false)
{
  int N=final_vectors.length();
  if (N<2) PLERROR("PowerIteration: final_vectors.length_ = %d should be >= 2",N);
  Vec previous=final_vectors(0);
  Vec current=final_vectors(1);
  previous << x0; 
  real max_x = max(previous);
  if (max_x<0) previous*=-1;
  real current_Rq, previous_Rq=0;
  max_x = max(current);
  if (max_x<0) current*=-1;
  normalize(previous);
  int it=1;
  Vec Ax = x0; 
  real current_error =0;
  for (;it<=n_iterations;it++)
  {
    int CGniter = max_n_cg_iter;
    real residue = RTolerance;
    current << previous;
    bool success=SolveLinearSymmSystemByCG(A, current, previous, CGniter, residue, regularizer);
    if (verbose)
    {
      if (success)
        cout << "done CG in " << CGniter << " iterations with residue = " << residue << endl;
      else
        cout << "done incomplete CG in " << CGniter << " iterations with residue = " << residue << endl;
    }
    max_x = max(current);
    if (max_x<0) current*=-1;
    normalize(current);
    product(A, current,Ax);
    current_Rq = dot(current,Ax);
    current_error = norm(Ax);
    if (verbose)
      cout << "at iteration " << it << ", R(A,x) = " << current_Rq << ", |Ax|/|x| = " 
           << current_error << endl;
    if (verbose && current_Rq > (1+RTolerance)*previous_Rq)
      PLWARNING("InversePowerIteration: something strange, x'Ax/x'x is decreasing %g->%g",
              previous_Rq, current_Rq);
    if (it>=N && 1 - current_Rq / previous_Rq  < RTolerance)
    {
      n_iterations = it;
      final_offset = it%N;
      if (verbose)
        cout << "inverse power iteration finishes with |Ax|/|x| = " << current_error << endl;
      return current_error;
    }
    previous_Rq = current_Rq;
    previous = current;
    current = final_vectors((it+1)%N);
  }
  final_offset = it%N;
  if (verbose)
    cout << "power iteration finishes FAILING with |Ax|/|x| = " << current_error << endl;
  return current_error;
}

template<class MatT>
real findSmallestEigenPairOfSymmMat(MatT& A, Vec x, 
                                    real error_tolerance=1e-3,
                                    real improvement_tolerance=1e-4, 
                                    int max_n_cg_iter=0, 
                                    int max_n_power_iter=0, bool verbose=false)
{
  int n=A.length();
  int n_try=5;

  if (max_n_cg_iter==0)
    max_n_cg_iter = 5+int(pow(double(n),0.3));
  if (max_n_power_iter==0)
    max_n_power_iter = 5+int(log(n));

  Mat try_solutions(n_try,n);
  Mat kernel_solutions = x.toMat(1,n);
  Vec Ax(n);

  int max_iter = int(sqrt(max_n_power_iter));
  real err=1e30, prev_err=1e30;
  int nrepeat=0;
  do {
    int n_iter=max_iter;
    int offs=0;
    real l0 = InversePowerIteration(A,x,n_iter,max_n_cg_iter,
                                    improvement_tolerance,
                                    try_solutions,offs,error_tolerance,verbose);
    if (verbose)
      cout << "got smallest eigenvalue = " << l0
           << " in " << n_iter << " iterations" << endl;

    if (verbose)
    {
      n_try = try_solutions.length();
      for (int i=0;i<n_try;i++)
        {
          cout << "for solution " << i << endl;
          Vec y = try_solutions(i);
          normalize(y);
          product(A, y,Ax);
          prev_err=err;
          err = norm(Ax);
          cout << "|A y| = " << err << endl;
          cout << "y. A y = " << dot(y,Ax) << endl;
        }
    }

    Vec evalues(n_try);
    Mat evectors(n_try,n_try);
    diagonalizeSubspace(A, try_solutions, Ax, kernel_solutions, evalues, evectors);
    product(A, x,Ax);
    err = norm(Ax);
    if (verbose)
      cout << "after diagonalization, err = " << err << endl;
    nrepeat++;
  } while (err > error_tolerance && n_try>=2 && 
           prev_err/err-1>improvement_tolerance && nrepeat<max_iter);
  if (verbose)
    cout << "return from findSmallestEigenPairOfSymmMat with err=" << err << endl;
}


template<class MatT>
int SymmMatNullSpaceByInversePowerIteration(MatT &A, Mat solutions,
                                            Vec normsAx, Vec xAx,  
                                            real error_tolerance=1e-3,
                                            real improvement_tolerance=1e-4, 
                                            int max_n_cg_iter=0, 
                                            int max_n_power_iter=0, bool verbose=false)

{  
  int n=A.length();
  int n_soln=normsAx.length();
  solutions.resize(n_soln,n);
  if (max_n_cg_iter==0)
    max_n_cg_iter = 5+int(pow(double(n),0.3));
  if (max_n_power_iter==0)
    max_n_power_iter = 5+int(log(n));
  Vec Ax(n);

  MatTPlusSumSquaredVec<MatT> B(A,n_soln);
  Vec sy(n);
  fill_random_uniform(sy);
  Mat largest_evecs(3,n);
  int offs;
  int n_iter = MIN(max_n_power_iter,20);
  real largest_evalue = PowerIteration(A, sy, n_iter ,1e-3,
                                       largest_evecs, offs,verbose);
  if (verbose)
    cout << "largest evalue(B) = " << largest_evalue << endl;
  for (int i=0;i<n_soln;i++)
  {
    Vec y = solutions(i);
    if (i==0)
      y.fill(1);
    else
      fill_random_uniform(y,0.1,0.5);
    real residue=findSmallestEigenPairOfSymmMat(B, y, 
                                                error_tolerance,
                                                improvement_tolerance,
                                                max_n_cg_iter,
                                                max_n_power_iter,verbose);
    if (verbose)
    {
      cout << "found vector with |B y| = " << residue << endl;
      cout << "|y| = " << norm(y) << endl;
      product(A, y,Ax);
      cout << "****** |A y| = " << norm(Ax) << endl;
    }
    multiply(y,largest_evalue,sy);
    B.squaredVecAcc(sy);
  }

  int n_s=GramSchmidtOrthogonalization(solutions);
  solutions = solutions.subMatRows(0,n_s);
  if (n_s<n_soln && verbose)
    cout << "found only " << n_s << " independent solutions out of " 
         << n_soln << " requested" << endl;
  for (int i=0;i<n_s;i++)
  {
    Vec xi = solutions(i);
    product(A, xi,Ax);
    real err = dot(xi,Ax);
    real normAx = norm(Ax);
    normsAx[i]=normAx;
    xAx[i]=err;
    if (verbose)
      cout << "for " << i << "-th solution x: x'Ax = " << err << ", |Ax|/|x|= "
           << normAx << endl;
  }
  return n_s;
}

template<class MatT>
int GDFindSmallEigenPairs(MatT& A,Mat X,
                           bool diagonalize_in_the_end=true,
                           real tolerance=1e-6,
                           int n_epochs=0,
                           real learning_rate=0,
                           int normalize_every=0,
                           real decrease_factor=0,
                           bool verbose=false)
{
  int n=A.length();
  int m=X.length();
  if (n_epochs==0)
    n_epochs = 1000+10*int(sqrt(n));
  Vec Ax(n);
  real err_tolerance, sum_norms, actual_err;
  if (learning_rate==0)
  {
    fill_random_uniform(Ax,-1,1);
    Mat large_vectors(3,n);
    int offs;
    int n_iter = 10+int(sqrt(log(double(n))));
    real max_eigen_value = 
      PowerIteration(A, Ax, n_iter,1e-3,large_vectors,offs,verbose);
    learning_rate = 2.0/max_eigen_value;
    if (verbose)
    {
      cout << "setting initial learning rate = 2/max_eigen_value = 2/" 
           << max_eigen_value << " = " << learning_rate << endl;
    }
    err_tolerance = tolerance * max_eigen_value;
  }
  else err_tolerance = tolerance;
  real prev_err=1e30;
  int it=0;
  for (;it<n_epochs;it++)
  {
    real learning_rate = learning_rate / (1+it*decrease_factor);
    real err=0;
    for (int i=0;i<n;i++)
    {
      real* xi = &X(0,i);
      for (int d=0;d<m;d++, xi+=n)
      {
        real gradient = matRowDotVec(A, i,X(d));
        *xi -= learning_rate * gradient;
        err += gradient * *xi;
      }
    }
    if (verbose)
    {
      cout << "at iteration " << it << " of gradient descent, est. err.= " << err << endl;
      cout << "lrate = " << learning_rate << endl;
    }
    if (tolerance>0)
    {
      sum_norms = 0;
      for (int d=0;d<m;d++)
        sum_norms += norm(X(d));
      actual_err = err / (m*sum_norms);
      if (actual_err<err_tolerance) break;
    }
    if (err>prev_err)
      cout << "at iteration " << it << " of gradient descent, est. err.= " << err 
           << " > previous error = " << prev_err << " ; lrate = " << learning_rate << endl;
    if (verbose)
      for (int d=0;d<m;d++)
        cout << "norm(x[" << d << "])=" << norm(X(d)) << endl;
    if (normalize_every!=0 && it%normalize_every==0)
    {
      int new_m=GramSchmidtOrthogonalization(X,1e-9);
      for (int e=new_m;e<m;e++)
        fill_random_uniform(X(e)); 
      if (verbose)
      {
        real C=0;
        for (int d=0;d<m;d++)
          {
            Vec xd = X(d);
            product(A, xd,Ax);
            C += dot(xd,Ax);
            cout << "for " << d << ", |Ax|/|x| = " << norm(Ax) << endl;
          }
        cout << "actual cost = " << 0.5*C << endl;
      }
    }
  }
  if (diagonalize_in_the_end)
  {
    Mat diagonalized_solutions(m,n);
    Mat subspace_evectors(m,m);
    Vec subspace_evalues(m);
    diagonalizeSubspace(A,X,Ax,diagonalized_solutions,subspace_evalues,subspace_evectors);
    X << diagonalized_solutions;
  }
  return it;
}

template<class MatT>
int metricMultiDimensionalScaling(MatT& square_distances,Mat embedding, int max_n_eigen_iter=300)
{
  int n=square_distances.length();
  long int m=embedding.width();
  if (embedding.length()!=n)
    PLERROR("MetricMultiDimensionalScaling: expected embedding.length()==square_distances.length(), got %d!=%d",
          embedding.length(),n);
  if (square_distances.width()!=n)
    PLERROR("MetricMultiDimensionalScaling: expected square_distances a square matrix, got %d x %d",
          n,square_distances.width());
  if (square_distances.size()!=n*n)
    PLWARNING("MetricMultiDimensionalScaling: current implementation does not seem to work well for sparse distance matrices!");  
  Vec avg_across_rows(n), avg_across_columns(n);
  averageAcrossRowsAndColumns(square_distances, avg_across_rows, avg_across_columns, true);
  real overall_avg=mean(avg_across_rows);
  avg_across_rows *= -1; 
  avg_across_columns *= -1; 
  addToRows(square_distances, avg_across_rows,true); 
  addToColumns(square_distances, avg_across_rows,true); 
  square_distances.add(overall_avg,true);
  square_distances*=-0.5;

  Vec e_values(m);
  Mat e_vectors(m,n);
  int err=eigenSparseSymmMat(square_distances, e_values, 
                             e_vectors, m, max_n_eigen_iter);
  if (err!=0)
    return err;
  for (int j=0;j<m;j++)
  {
    real eval_j = e_values[j];
    if (eval_j<0)
      PLERROR("metricMultiDimensionalScaling::the matrix of dot-products is not positive-definite!, evalue=%g",eval_j);
    real scale = sqrt(eval_j);
    Vec feature_j = e_vectors(j);
    feature_j *= scale;
    embedding.column(j) << feature_j;
  }
  return 0;
}


} // end of namespace PLearn

#endif

---