Package org.apache.mahout.math.decomposer.lanczos

Source Code of org.apache.mahout.math.decomposer.lanczos.LanczosSolver$Scale

/**
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements.  See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License.  You may obtain a copy of the License at
*
*     http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/

package org.apache.mahout.math.decomposer.lanczos;


import org.apache.mahout.math.Matrix;
import org.apache.mahout.math.Vector;
import org.apache.mahout.math.VectorIterable;
import org.apache.mahout.math.function.DoubleFunction;
import org.apache.mahout.math.function.PlusMult;
import org.apache.mahout.math.matrix.DoubleMatrix1D;
import org.apache.mahout.math.matrix.DoubleMatrix2D;
import org.apache.mahout.math.matrix.linalg.EigenvalueDecomposition;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import java.util.EnumMap;
import java.util.Map;

/**
* <p>Simple implementation of the <a href="http://en.wikipedia.org/wiki/Lanczos_algorithm">Lanczos algorithm</a> for
* finding eigenvalues of a symmetric matrix, applied to non-symmetric matrices by applying Matrix.timesSquared(vector)
* as the "matrix-multiplication" method.</p>
* <p>
* To avoid floating point overflow problems which arise in power-methods like Lanczos, an initial pass is made
* through the input matrix to
* <ul>
*   <li>generate a good starting seed vector by summing all the rows of the input matrix, and</li>
*   <li>compute the trace(inputMatrix<sup>t</sup>*matrix)
* </ul>
* </p>
* <p>
* This latter value, being the sum of all of the singular values, is used to rescale the entire matrix, effectively
* forcing the largest singular value to be strictly less than one, and transforming floating point <em>overflow</em>
* problems into floating point <em>underflow</em> (ie, very small singular values will become invisible, as they
* will appear to be zero and the algorithm will terminate).
* </p>
* <p>This implementation uses {@link org.apache.mahout.math.matrix.linalg.EigenvalueDecomposition} to do the
* eigenvalue extraction from the small (desiredRank x desiredRank) tridiagonal matrix.  Numerical stability is
* achieved via brute-force: re-orthogonalization against all previous eigenvectors is computed after every pass.
* This can be made smarter if (when!) this proves to be a major bottleneck.  Of course, this step can be parallelized
* as well.
* </p>
*/
public class LanczosSolver {

  private static final Logger log = LoggerFactory.getLogger(LanczosSolver.class);

  public static final double SAFE_MAX = 1.0e150;

  public enum TimingSection {
    ITERATE, ORTHOGANLIZE, TRIDIAG_DECOMP, FINAL_EIGEN_CREATE
  }

  private final Map<TimingSection, Long> startTimes = new EnumMap<TimingSection, Long>(TimingSection.class);
  private final Map<TimingSection, Long> times = new EnumMap<TimingSection, Long>(TimingSection.class);

  private static final class Scale implements DoubleFunction {
    private final double d;

    private Scale(double d) {
      this.d = d;
    }

    @Override
    public double apply(double arg1) {
      return arg1 * d;
    }
  }

  public void solve(LanczosState state,
                    int desiredRank) {
    solve(state, desiredRank, false);
  }

  public void solve(LanczosState state,
                    int desiredRank,
                    boolean isSymmetric) {
    VectorIterable corpus = state.getCorpus();
    log.info("Finding {} singular vectors of matrix with {} rows, via Lanczos",
        desiredRank, corpus.numRows());
    int i = state.getIterationNumber();
    Vector currentVector = state.getBasisVector(i - 1);
    Vector previousVector = state.getBasisVector(i - 2);
    double beta = 0;
    Matrix triDiag = state.getDiagonalMatrix();
    while (i < desiredRank) {
      startTime(TimingSection.ITERATE);
      Vector nextVector = isSymmetric ? corpus.times(currentVector) : corpus.timesSquared(currentVector);
      log.info("{} passes through the corpus so far...", i);
      if(state.getScaleFactor() <= 0) {
        state.setScaleFactor(calculateScaleFactor(nextVector));
      }
      nextVector.assign(new Scale(1.0 / state.getScaleFactor()));
      if(previousVector != null) {
        nextVector.assign(previousVector, new PlusMult(-beta));
      }
      // now orthogonalize
      double alpha = currentVector.dot(nextVector);
      nextVector.assign(currentVector, new PlusMult(-alpha));
      endTime(TimingSection.ITERATE);
      startTime(TimingSection.ORTHOGANLIZE);
      orthoganalizeAgainstAllButLast(nextVector, state);
      endTime(TimingSection.ORTHOGANLIZE);
      // and normalize
      beta = nextVector.norm(2);
      if (outOfRange(beta) || outOfRange(alpha)) {
        log.warn("Lanczos parameters out of range: alpha = {}, beta = {}.  Bailing out early!",
            alpha, beta);
        break;
      }
      nextVector.assign(new Scale(1 / beta));
      state.setBasisVector(i, nextVector);
      previousVector = currentVector;
      currentVector = nextVector;
      // save the projections and norms!
      triDiag.set(i - 1, i - 1, alpha);
      if (i < desiredRank - 1) {
        triDiag.set(i - 1, i, beta);
        triDiag.set(i, i - 1, beta);
      }
      state.setIterationNumber(++i);
    }
    startTime(TimingSection.TRIDIAG_DECOMP);

    log.info("Lanczos iteration complete - now to diagonalize the tri-diagonal auxiliary matrix.");
    // at this point, have tridiag all filled out, and basis is all filled out, and orthonormalized
    EigenvalueDecomposition decomp = new EigenvalueDecomposition(triDiag);

    DoubleMatrix2D eigenVects = decomp.getV();
    DoubleMatrix1D eigenVals = decomp.getRealEigenvalues();
    endTime(TimingSection.TRIDIAG_DECOMP);
    startTime(TimingSection.FINAL_EIGEN_CREATE);
    for (int row = 0; row < i; row++) {
      Vector realEigen = null;
      // the eigenvectors live as columns of V, in reverse order.  Weird but true.
      DoubleMatrix1D ejCol = eigenVects.viewColumn(i - row - 1);
      int size = Math.min(ejCol.size(), state.getBasisSize());
      for (int j = 0; j < size; j++) {
        double d = ejCol.get(j);
        Vector rowJ = state.getBasisVector(j);
        if(realEigen == null) {
          realEigen = rowJ.like();
        }
        realEigen.assign(rowJ, new PlusMult(d));
      }
      realEigen = realEigen.normalize();
      state.setRightSingularVector(row, realEigen);
      double e = eigenVals.get(row) * state.getScaleFactor();
      if(!isSymmetric) {
        e = Math.sqrt(e);
      }
      log.info("Eigenvector {} found with eigenvalue {}", row, e);
      state.setSingularValue(row, e);
    }
    log.info("LanczosSolver finished.");
    endTime(TimingSection.FINAL_EIGEN_CREATE);
  }

  protected double calculateScaleFactor(Vector nextVector) {
    return nextVector.norm(2);
  }

  private static boolean outOfRange(double d) {
    return Double.isNaN(d) || d > SAFE_MAX || -d > SAFE_MAX;
  }

  protected void orthoganalizeAgainstAllButLast(Vector nextVector, LanczosState state) {
    for (int i = 0; i < state.getIterationNumber(); i++) {
      Vector basisVector = state.getBasisVector(i);
      double alpha;
      if (basisVector == null || (alpha = nextVector.dot(basisVector)) == 0.0) {
        continue;
      }
      nextVector.assign(basisVector, new PlusMult(-alpha));
    }
  }

  private void startTime(TimingSection section) {
    startTimes.put(section, System.nanoTime());
  }

  private void endTime(TimingSection section) {
    if (!times.containsKey(section)) {
      times.put(section, 0L);
    }
    times.put(section, times.get(section) + System.nanoTime() - startTimes.get(section));
  }

}
TOP

Related Classes of org.apache.mahout.math.decomposer.lanczos.LanczosSolver$Scale

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.