0
\$\begingroup\$

My C++20 N-dimensional matrix project now supports basic linear algebra operations: https://github.com/frozenca/Ndim-Matrix

Today I want to get some reviews on computing eigenvalues and eigenvectors. In particular, I don't like so much duplicate codes. I think other parts are okay.. (hope?)

There are two functions: Computes eigenvalues only, and computes both eigenvalues and corresponding eigenvectors. These two functions share so much code, I want to know how to fix nicely..

Full code: https://github.com/frozenca/Ndim-Matrix/blob/main/LinalgOps.h#L956

constexpr float tolerance_soft = 1e-6;
constexpr float tolerance_hard = 1e-10;
constexpr std::size_t max_iter = 100;
constexpr std::size_t local_iter = 15;

template <isScalar T>
decltype(auto) getSign (const T& t) {
    if constexpr (isComplex<T>) {
        if (std::abs(t) < tolerance_soft) {
            return static_cast<T>(1.0f);
        }
        return t / std::abs(t);
    } else {
        return std::signbit(t) ? -1.0f : +1.0f;
    }
};

template <typename Derived, isScalar T>
Vec<T> normalize(const MatrixBase<Derived, T, 1>& vec) {
    Vec<T> u = vec;
    auto u_norm = norm(u);
    if (u_norm < tolerance_hard) {
        return u;
    } else {
        u /= u_norm;
        return u;
    }
}

template <typename Derived, isScalar T>
bool isSubdiagonalNeglegible(MatrixBase<Derived, T, 2>& M, std::size_t idx) {
    if (std::abs(M[{idx + 1, idx}]) <
        tolerance_soft * (std::abs(M[{idx, idx}]) + std::abs(M[{idx + 1, idx + 1}]))) {
        M[{idx + 1, idx}] = T{0};
        return true;
    }
    return false;
}

template <isScalar T>
T computeShift(const T& a, const T& b, const T& c, const T& d) {
    auto tr = a + d;
    auto det = a * d - b * c;
    auto disc = std::sqrt(tr * tr - 4.0f * det);
    auto root1 = (tr + disc) / 2.0f;
    auto root2 = (tr - disc) / 2.0f;
    if (std::abs(root1 - d) < std::abs(root2 - d)) {
        return root1;
    } else {
        return root2;
    }
}

template <isScalar T, isReal U = RealTypeT<T>>
std::tuple<T, T, U> givensRotation(const T& a, const T& b) {
    if (b == T{0}) {
        return {getSign(a), 0, std::abs(a)};
    } else if (a == T{0}) {
        return {0, getSign(b), std::abs(b)};
    } else if (std::abs(a) > std::abs(b)) {
        auto t = b / a;
        auto u = getSign(a) * std::sqrt(1.0f + t * t);
        return {1.0f / u, t / u, std::abs(a * u)};
    } else {
        auto t = a / b;
        auto u = getSign(b) * std::sqrt(1.0f + t * t);
        return {t / u, 1.0f / u, std::abs(b * u)};
    }
}


// QR algorithm used in eigendecomposition.
// not to be confused with QR decomposition
template <typename Derived, isScalar U, isScalar T = CmpTypeT<U>> requires CmpTypeTo<U, T>
std::vector<T> QRIteration(const MatrixBase<Derived, U, 2>& mat) {
    std::size_t iter = 0;
    std::size_t total_iter = 0;
    std::size_t n = mat.dims(0);

    auto conjif = [&](const auto& v) {
        if constexpr (isComplex<U>) {
            return conj(v);
        } else {
            return v;
        }
    };

    Mat<T> M = mat;
    std::size_t p = n - 1;
    while (true) {
        while (p > 0) {
            if (!isSubdiagonalNeglegible(M, p - 1)) {
                break;
            }
            iter = 0;
            --p;
        }
        if (p == 0) {
            break;
        }
        if (++iter > local_iter) {
            break;
        }
        if (++total_iter > max_iter) {
            break;
        }
        std::size_t top = p - 1;
        while (top > 0 && !isSubdiagonalNeglegible(M, top - 1)) {
            --top;
        }
        auto shift = computeShift(M[{p - 1, p - 1}], M[{p - 1, p}],
                                  M[{p, p - 1}], M[{p, p}]);

        // initial Givens rotation
        auto x = M[{top, top}] - shift;
        auto y = M[{top + 1, top}];
        auto [c, s, r] = givensRotation(x, y);
        Mat<T> R {{c, -s},
                  {s, c}};
        Mat<T> RT{{c, s},
                  {-s, c}};
        if (r > tolerance_hard) {
            auto Sub1 = M.submatrix({top, top}, {top + 2, n});
            Sub1 = dot(conjif(RT), Sub1);
            std::size_t bottom = std::min(top + 3, p + 1);
            auto Sub2 = M.submatrix({0, top}, {bottom, top + 2});
            Sub2 = dot(Sub2, R);
        }
        for (std::size_t k = top + 1; k < p; ++k) {
            x = M[{k, k - 1}];
            y = M[{k + 1, k - 1}];
            std::tie(c, s, r) = givensRotation(x, y);
            if (r > tolerance_hard) {
                M[{k, k - 1}] = r;
                M[{k + 1, k - 1}] = T{0};
                R[{0, 0}] = RT[{0, 0}] = R[{1, 1}] = RT[{1, 1}] = c;
                R[{0, 1}] = RT[{1, 0}] = -s;
                R[{1, 0}] = RT[{0, 1}] = s;

                auto Sub1 = M.submatrix({k, k}, {k + 2, n});
                Sub1 = dot(conjif(RT), Sub1);
                std::size_t bottom = std::min(k + 3, p + 1);
                auto Sub2 = M.submatrix({0, k}, {bottom, k + 2});
                Sub2 = dot(Sub2, R);
            }
        }

    }
    std::vector<T> res;
    for (std::size_t k = 0; k < n; ++k) {
        res.push_back(M[{k, k}]);
    }
    return res;
}

template <typename Derived, typename Derived2, isScalar T>
Mat<T> computeEigenvectors(const MatrixBase<Derived, T, 2>& M,
                           const MatrixBase<Derived2, T, 2>& Q) {
    std::size_t n = M.dims(0);
    Mat<T> X = identity<T>(n);

    for (std::size_t k = n - 1; k < n; --k) {
        for (std::size_t i = k - 1; i < n; --i) {
            X[{i, k}] -= M[{i, k}];
            if (k - i > 1 && k - i - 1 < n) {
                auto row_vec = M.row(i).submatrix(i + 1, k);
                auto col_vec = X.col(k).submatrix(i + 1, k);
                X[{i, k}] -= dot(row_vec, col_vec);
            }
            auto z = M[{i, i}] - M[{k, k}];
            if (z == T{0}) {
                z = static_cast<T>(tolerance_hard);
            }
            X[{i, k}] /= z;
        }
    }
    return dot(Q, X);
}

// QR algorithm used in eigendecomposition.
// not to be confused with QR decomposition
template <typename Derived, typename Derived2,
        isScalar U, isScalar T = CmpTypeT<U>> requires CmpTypeTo<U, T>
std::vector<std::pair<T, Vec<T>>> QRIterationWithVec(const MatrixBase<Derived, U, 2>& mat,
                                  const MatrixBase<Derived2, U, 2>& V) {
    std::size_t iter = 0;
    std::size_t total_iter = 0;
    std::size_t n = mat.dims(0);

    auto conjif = [&](const auto& v) {
        if constexpr (isComplex<U>) {
            return conj(v);
        } else {
            return v;
        }
    };

    Mat<T> M = mat;
    std::size_t p = n - 1;
    Mat<T> Q = V;
    while (true) {
        while (p > 0) {
            if (!isSubdiagonalNeglegible(M, p - 1)) {
                break;
            }
            iter = 0;
            --p;
        }
        if (p == 0) {
            break;
        }
        if (++iter > 20) {
            break;
        }
        if (++total_iter > max_iter) {
            break;
        }
        std::size_t top = p - 1;
        while (top > 0 && !isSubdiagonalNeglegible(M, top - 1)) {
            --top;
        }
        auto shift = computeShift(M[{p - 1, p - 1}], M[{p - 1, p}],
                                  M[{p, p - 1}], M[{p, p}]);

        // initial Givens rotation
        auto x = M[{top, top}] - shift;
        auto y = M[{top + 1, top}];
        auto [c, s, r] = givensRotation(x, y);
        Mat<T> R {{c, -s},
                  {s, c}};
        Mat<T> RT{{c, s},
                  {-s, c}};
        if (r > tolerance_hard) {
            auto Sub1 = M.submatrix({top, top}, {top + 2, n});
            Sub1 = dot(conjif(RT), Sub1);
            std::size_t bottom = std::min(top + 3, p + 1);
            auto Sub2 = M.submatrix({0, top}, {bottom, top + 2});
            Sub2 = dot(Sub2, R);
            auto QSub2 = Q.submatrix({0, top}, {n, top + 2});
            QSub2 = dot(QSub2, R);
        }
        for (std::size_t k = top + 1; k < p; ++k) {
            x = M[{k, k - 1}];
            y = M[{k + 1, k - 1}];
            std::tie(c, s, r) = givensRotation(x, y);
            if (r > tolerance_hard) {
                M[{k, k - 1}] = r;
                M[{k + 1, k - 1}] = T{0};
                R[{0, 0}] = RT[{0, 0}] = R[{1, 1}] = RT[{1, 1}] = c;
                R[{0, 1}] = RT[{1, 0}] = -s;
                R[{1, 0}] = RT[{0, 1}] = s;

                auto Sub1 = M.submatrix({k, k}, {k + 2, n});
                Sub1 = dot(conjif(RT), Sub1);
                std::size_t bottom = std::min(k + 3, p + 1);
                auto Sub2 = M.submatrix({0, k}, {bottom, k + 2});
                Sub2 = dot(Sub2, R);
                auto QSub2 = Q.submatrix({0, k}, {n, k + 2});
                QSub2 = dot(QSub2, R);
            }
        }
    }
    std::vector<std::pair<T, Vec<T>>> res;
    auto eV = computeEigenvectors(M, Q);
    for (std::size_t k = 0; k < n; ++k) {
        res.emplace_back(M[{k, k}], normalize(eV.col(k)));
    }
    return res;
}

} // anonymous namespace

template <typename Derived, isScalar U, isScalar T = CmpTypeT<U>> requires CmpTypeTo<U, T>
std::vector<T> eigenval(const MatrixBase<Derived, U, 2>& M) {
    std::size_t n = M.dims(0);
    std::size_t C = M.dims(1);
    if (n != C) {
        throw std::invalid_argument("Not a square Matrix, cannot compute eigenvalues");
    }

    if (n == 1) { // 1 x 1
        return {M[{0, 0}]};
    } else if (n == 2) { // 2 x 2
        return eigenTwo(M);
    } else if (n == 3) { // 3 x 3
        return eigenThree(M);
    } else { // for 4 x 4 we need advanced algorithm
        auto H = Hessenberg(M);
        return QRIteration(H);
    }
}

template <typename Derived, isScalar U, isScalar T = CmpTypeT<U>> requires CmpTypeTo<U, T>
std::vector<std::pair<T, Vec<T>>> eigenvec(const MatrixBase<Derived, U, 2>& M) {
    std::size_t n = M.dims(0);
    std::size_t C = M.dims(1);
    if (n != C) {
        throw std::invalid_argument("Not a square Matrix, cannot compute eigenvalues");
    }

    if (n == 1) { // 1 x 1
        auto val = M[{0, 0}];
        Vec<T> vec {T{1}};
        return {{val, vec}};
    } else if (n == 2) { // 2 x 2
        return eigenVecTwo(M);
    } else if (n == 3) { // 3 x 3
        return eigenVecThree(M);
    } else { // for 4 x 4 we need advanced algorithm
        auto [H, V] = HessenbergWithVec(M);
        return QRIterationWithVec(H, V);
    }
}

```
\$\endgroup\$
2
  • \$\begingroup\$ Can you edit your question to include the headers necessary for this to compile? \$\endgroup\$ – 1201ProgramAlarm May 8 at 2:28
  • \$\begingroup\$ @1201ProgramAlarm Very sorry, there are too many dependencies... to compile, you should git clone https://github.com/frozenca/Ndim-Matrix.git and add #include "Matrix.h". The test.cpp shows the usage of the code I posted. \$\endgroup\$ – frozenca May 8 at 2:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.