Уменьшенное разложение QR с собственным - PullRequest
Новая
       4

Уменьшенное разложение QR с собственным

0 голосов
/ 07 января 2020

Я пытаюсь решить классическую проблему наименьших квадратов argmin (Ax - b)^2, используя Eigen. Я знаю, что система переопределена (то есть m >= n, где A равно m x n), и что A имеет полный ранг (n). Я решил использовать разреженную QR-часть Eigen, поскольку A сам по себе является разреженным. Я дошел до этого: 

#include <Eigen/Sparse>
#include <Eigen/SparseQR>

void solve_least_squares(const Eigen::SparseMatrix<double>& matrix,
                         const Eigen::VectorXd& rhs)
{
  int m = matrix.rows();
  int n = matrix.cols();

  assert(m >= n);
  assert(matrix.isCompressed());

  Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> decomposition(matrix);

  Eigen::VectorXd t = (decomposition.matrixQ().transpose() * rhs).topRows(n);

  std::cout << t << std::endl;

  Eigen::VectorXd s = decomposition.matrixR().topRows(n).triangularView<Eigen::Upper>().solve(t);

  return decomposition.colsPermutation().inverse()*s;
}

Но мне интересно, является ли это наиболее эффективной реализацией: во-первых, Эйген, похоже, определяет полное QR-разложение, а не сокращенное ( Q это m x m, а не m x n). Это кажется расточительным, поскольку мне нужны только верхние n строки.

В плотном случае есть любой пример в классе HouseholderQR : 

MatrixXf A(MatrixXf::Random(5,3)), thinQ(MatrixXf::Identity(5,3)), Q;
A.setRandom();
HouseholderQR<MatrixXf> qr(A);
Q = qr.householderQ();
thinQ = qr.householderQ() * thinQ; // <- here
std::cout << "The complete unitary matrix Q is:\n" << Q << "\n\n";
std::cout << "The thin matrix Q is:\n" << thinQ << "\n\n";

Здесь умножение матриц используется для получения уменьшенного Q, что кажется даже более неэффективным, чем нарезка полной матрицы.

Поскольку разложение (собственного) SVD (плотного) SVD дает возможность вычислять тонкий SVD, есть ли подобная опция для разложения QR?

1 Ответ

1 голос
/ 07 января 2020

Внутри разложения коэффициент Q сохраняется как последовательность векторов Домохозяина, а метод matrixQ() по существу возвращает ссылку на эту матрицу (которая перегружает умножение таким образом, как если бы оно умножалось на фактическую матрицу). , При хранении в качестве матрицы домохозяев, не имеет значения, если Q представляет полную или тонкую часть Q (на самом деле, умножение вектора на полную матрицу может быть эффективно выполнено на месте).

Если Если вы хотите решить свою линейную систему, вы должны просто написать 

Eigen::VectorXd s = decomposition.solve(rhs);

Кстати: если вам хорошо с потерей какой-то точности (и ваша матрица достаточно хорошо подготовлена), то вы (в большинстве случаев) ) быстрее, решая нормальное уравнение (запись Matlab): 

x = (A'*A) \ (A'*b)

и используя разреженное разложение Холецкого (например, Eigen::SimplicialLLT или Eigen::SimplicialLDLT). Но сравните это с вашими фактическими данными и проверьте, достаточна ли точность для вашего варианта использования (возможно, после этапа уточнения, в котором повторно используется разложение Холецкого).

...