max_eigenvector.h

Go to the documentation of this file.

/*
    Copyright (C) 2005 - 2007-10-08, Hammersmith Imanet Ltd
    This file is part of STIR.

    SPDX-License-Identifier: Apache-2.0

    See STIR/LICENSE.txt for details
*/
#ifndef __stir_numerics_max_eigenvector_H__
#define __stir_numerics_max_eigenvector_H__

#include "stir/numerics/MatrixFunction.h"
#include "stir/numerics/norm.h"
#include "stir/more_algorithms.h"
#include "stir/Succeeded.h"
#include "stir/error.h"

START_NAMESPACE_STIR

template <class elemT>
inline Succeeded
absolute_max_eigenvector_using_power_method(elemT& max_eigenvalue,
                                            Array<1, elemT>& max_eigenvector,
                                            const Array<2, elemT>& m,
                                            const Array<1, elemT>& start,
                                            const double tolerance = .01,
                                            const unsigned long max_num_iterations = 10000UL)
{
  assert(m.is_regular());
  if (m.size() == 0)
    {
      max_eigenvalue = 0;
      max_eigenvector = start;
      return Succeeded::yes;
    }

  const double tolerance_squared = square(tolerance);
  Array<1, elemT> current = start;
  current /= (*abs_max_element(current.begin(), current.end()));
  unsigned long remaining_num_iterations = max_num_iterations;

  double change;
  do
    {
      max_eigenvector = matrix_multiply(m, current);
      const elemT norm_factor = *abs_max_element(max_eigenvector.begin(), max_eigenvector.end());
      max_eigenvector /= norm_factor;
      change = norm_squared(max_eigenvector - current);
      current = max_eigenvector;
      --remaining_num_iterations;
  } while (change > tolerance_squared && remaining_num_iterations != 0);

  current /= static_cast<elemT>(norm(current));
  max_eigenvector = matrix_multiply(m, current);
  // compute eigenvalue using Rayleigh quotient
  max_eigenvalue = static_cast<elemT>(inner_product(current, max_eigenvector) / norm_squared(current));
  max_eigenvector /= static_cast<elemT>(norm(max_eigenvector));

  return remaining_num_iterations == 0 ? Succeeded::no : Succeeded::yes;

  /*norm( max_eigenvector*max_eigenvalue -
               matrix_multiply(m, max_eigenvector));*/
}

template <class elemT>
inline Succeeded
absolute_max_eigenvector_using_shifted_power_method(elemT& max_eigenvalue,
                                                    Array<1, elemT>& max_eigenvector,
                                                    const Array<2, elemT>& m,
                                                    const Array<1, elemT>& start,
                                                    const elemT shift,
                                                    const double tolerance = .03,
                                                    const unsigned long max_num_iterations = 10000UL)
{
  if (m.get_min_index() != 0 || m[0].get_min_index() != 0)
    error("absolute_max_eigenvector_using_shifted_power_method:\n"
          "  implementation needs work for indices that don't start from 0. sorry");

  Succeeded success
      = absolute_max_eigenvector_using_power_method(max_eigenvalue,
                                                    max_eigenvector,
                                                    // sadly need to explicitly convert result of subtraction back to Array
                                                    Array<2, elemT>(m - diagonal_matrix(static_cast<unsigned>(m.size()), shift)),
                                                    start,
                                                    tolerance,
                                                    max_num_iterations);
  max_eigenvalue += shift;
  return success;
}

template <class elemT>
inline Succeeded
max_eigenvector_using_power_method(elemT& max_eigenvalue,
                                   Array<1, elemT>& max_eigenvector,
                                   const Array<2, elemT>& m,
                                   const Array<1, elemT>& start,
                                   const double tolerance = .03,
                                   const unsigned long max_num_iterations = 10000UL)
{

  Succeeded success
      = absolute_max_eigenvector_using_power_method(max_eigenvalue, max_eigenvector, m, start, tolerance, max_num_iterations);
  if (success == Succeeded::no)
    return Succeeded::no;

  if (max_eigenvalue >= 0) // TODO would need to take real value for complex case
    return Succeeded::yes;

  // we found a negative eigenvalue
  // try again with a shift equal to the previously found max_eigenvalue
  // this shift will effectively put that eigenvalue to 0 during the
  // power method iterations
  // also it will make all eigenvalues positive (as we subtract the
  // smallest negative eigenvalue)
  success = absolute_max_eigenvector_using_shifted_power_method(
      max_eigenvalue, max_eigenvector, m, start, max_eigenvalue, tolerance, max_num_iterations);
  if (success == Succeeded::no)
    return Succeeded::no;

  assert(max_eigenvalue >= 0);
  return Succeeded::yes;
}

END_NAMESPACE_STIR

#endif