NicheLibrary

This documentation is automatically generated by NotLeonian/competitive-verifier (forked from competitive-verifier/competitive-verifier)

View the Project on GitHub NotLeonian/NicheLibrary

:heavy_check_mark: 行列の階数と動的な 1 行・1 列更新 (math/matrix/dynamic-matrix-rank.hpp)

概要

  • 体上の $r \times c$ 行列を前処理し、現在の階数と $A + uv^\top$ の階数を求める。
  • さらに、1 行差し替え・1 列差し替え・外積 1 項更新を内部状態に反映できる。
  • 現在の行列は左右の階数分解と片側逆元で保持する。

使い方

  • DynamicMatrixRank()
    • 空に構築する。後で build(matrix) を呼ぶ。
  • DynamicMatrixRank(const std::vector<std::vector<T>>& matrix)
    • matrix で構築する。
    • 前提：matrix は長方形であり、T は除算ができる。
  • void build(const std::vector<std::vector<T>>& matrix)
    • matrix を現在の行列として前処理し直す。
  • void build()
    • 現在保持している行列から前処理し直す。
  • int rank() const
    • 現在の行列の階数を返す。
  • std::vector<T> get_row(int row_index) const
    • 現在の row_index 行目を返す。
    • 前提：$0\le \mathrm{row_index}<r$（$r$ は行数）。
  • std::vector<T> get_column(int column_index) const
    • 現在の column_index 列目を返す。
    • 前提：$0\le \mathrm{column_index}<c$（$c$ は列数）。
  • std::vector<std::vector<T>> materialize_matrix() const
    • 現在の行列を密行列として返す。
  • int rank_after_rank_one_update(const std::vector<T>& column_vector, const std::vector<T>& row_vector) const
    • $A+CR^{\top}$ の階数を返す（ただし、$C=\mathrm{column_vector},\,R=\mathrm{row_vector}$）。
    • 前提：column_vector の長さは行数、row_vector の長さは列数。
    • 備考：内部状態は変更しない。
  • int rank_after_row_replacement(int row_index, const std::vector<T>& new_row) const
    • row_index 行目を new_row に差し替えた行列の階数を返す。
    • 前提：$0\le \mathrm{row_index}<r$（$r$ は行数）、new_row の長さは列数。
    • 備考：内部状態は変更しない。
  • int rank_after_column_replacement(int column_index, const std::vector<T>& new_column) const
    • column_index 列目を new_column に差し替えた行列の階数を返す。
    • 前提：$0\le \mathrm{column_index}<c$（$c$ は列数）、new_column の長さは行数。
    • 備考：内部状態は変更しない。
  • int apply_rank_one_update(const std::vector<T>& column_vector, const std::vector<T>& row_vector)
    • $A+CR^{\top}$（ただし、$C=\mathrm{column_vector},\,R=\mathrm{row_vector}$）に内部状態を更新し、その階数を返す。
  • int apply_row_replacement(int row_index, const std::vector<T>& new_row)
    • row_index 行目を new_row に差し替え、変更後の階数を返す。
  • int apply_column_replacement(int column_index, const std::vector<T>& new_column)
    • column_index 列目を new_column に差し替え、変更後の階数を返す。

計算量

  • ここで $k$ は現在の階数である。
  • build：$O(rc\min(r, c) + k^2(r + c))$
  • rank：$O(1)$
  • get_row：$O(kc)$
  • get_column：$O(kr)$
  • materialize_matrix：$O(krc)$
  • rank_after_rank_one_update：$O(k(r + c))$
  • rank_after_row_replacement：$O(k(r + c))$
  • rank_after_column_replacement：$O(k(r + c))$
  • apply_rank_one_update：$O(k(r + c))$
  • apply_row_replacement：$O(k(r + c))$
  • apply_column_replacement：$O(k(r + c))$

Verified with

Code

#ifndef MATH_MATRIX_DYNAMIC_MATRIX_RANK_HPP
#define MATH_MATRIX_DYNAMIC_MATRIX_RANK_HPP

// 行列の階数を求め、外積 1 項更新後の階数を判定する。
// さらに、内部状態を O(k(r + c)) で更新しつつ変更後の階数を返せる。
// 現在の行列は左右の階数分解と片側逆元で保持する。
// 1 行差し替え・1 列差し替えはそれぞれ e_i (b-a_i)^T, (b-a_j) e_j^T に帰着する。
// 前処理は O(rc min(r, c) + k^2(r + c))、各更新・各判定は O(k(r + c)) である。
// T は体を成し、零判定と四則演算ができることを仮定する。

#include <cassert>
#include <vector>

template <class T> struct DynamicMatrixRank {
    int row_size = 0;
    int column_size = 0;
    int matrix_rank = 0;

    std::vector<std::vector<T>> column_space_basis;
    std::vector<std::vector<T>> row_space_basis;
    std::vector<std::vector<T>> column_space_left_inverse;
    std::vector<std::vector<T>> row_space_right_inverse;

    DynamicMatrixRank() = default;

    explicit DynamicMatrixRank(const std::vector<std::vector<T>> &matrix) {
        build(matrix);
    }

    void build(const std::vector<std::vector<T>> &matrix) {
        row_size = static_cast<int>(matrix.size());
        column_size = row_size == 0 ? 0 : static_cast<int>(matrix[0].size());
        for (int i = 1; i < row_size; ++i) {
            assert(static_cast<int>(matrix[i].size()) == column_size);
        }

        const std::vector<int> basis_columns = find_independent_columns(matrix);
        const std::vector<int> basis_rows =
            find_independent_columns(transpose_matrix(matrix));
        matrix_rank = static_cast<int>(basis_columns.size());
        assert(static_cast<int>(basis_rows.size()) == matrix_rank);

        std::vector<std::vector<T>> intersection_matrix(
            matrix_rank, std::vector<T>(matrix_rank, T()));
        for (int i = 0; i < matrix_rank; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                intersection_matrix[i][j] =
                    matrix[basis_rows[i]][basis_columns[j]];
            }
        }
        const std::vector<std::vector<T>> intersection_inverse =
            inverse_matrix(intersection_matrix);

        column_space_basis.assign(row_size, std::vector<T>(matrix_rank, T()));
        for (int i = 0; i < row_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                column_space_basis[i][j] = matrix[i][basis_columns[j]];
            }
        }

        row_space_basis.assign(matrix_rank, std::vector<T>(column_size, T()));
        for (int i = 0; i < matrix_rank; ++i) {
            for (int mid = 0; mid < matrix_rank; ++mid) {
                const T value = intersection_inverse[i][mid];
                if (value == T()) {
                    continue;
                }
                for (int j = 0; j < column_size; ++j) {
                    row_space_basis[i][j] += value * matrix[basis_rows[mid]][j];
                }
            }
        }

        column_space_left_inverse.assign(matrix_rank,
                                         std::vector<T>(row_size, T()));
        for (int i = 0; i < matrix_rank; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                column_space_left_inverse[i][basis_rows[j]] =
                    intersection_inverse[i][j];
            }
        }

        row_space_right_inverse.assign(column_size,
                                       std::vector<T>(matrix_rank, T()));
        for (int i = 0; i < matrix_rank; ++i) {
            row_space_right_inverse[basis_columns[i]][i] = T(1);
        }
    }

    void build() { build(materialize_matrix()); }

    int rank() const { return matrix_rank; }

    std::vector<T> get_row(int row_index) const {
        assert(0 <= row_index && row_index < row_size);
        std::vector<T> row(column_size, T());
        for (int i = 0; i < matrix_rank; ++i) {
            const T value = column_space_basis[row_index][i];
            if (value == T()) {
                continue;
            }
            for (int j = 0; j < column_size; ++j) {
                row[j] += value * row_space_basis[i][j];
            }
        }
        return row;
    }

    std::vector<T> get_column(int column_index) const {
        assert(0 <= column_index && column_index < column_size);
        std::vector<T> column(row_size, T());
        for (int i = 0; i < row_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                column[i] +=
                    column_space_basis[i][j] * row_space_basis[j][column_index];
            }
        }
        return column;
    }

    std::vector<std::vector<T>> materialize_matrix() const {
        std::vector<std::vector<T>> matrix(row_size,
                                           std::vector<T>(column_size, T()));
        for (int i = 0; i < row_size; ++i) {
            for (int mid = 0; mid < matrix_rank; ++mid) {
                const T value = column_space_basis[i][mid];
                if (value == T()) {
                    continue;
                }
                for (int j = 0; j < column_size; ++j) {
                    matrix[i][j] += value * row_space_basis[mid][j];
                }
            }
        }
        return matrix;
    }

    int rank_after_rank_one_update(const std::vector<T> &column_vector,
                                   const std::vector<T> &row_vector) const {
        return analyze_rank_one_update(column_vector, row_vector).next_rank;
    }

    int rank_after_row_replacement(int row_index,
                                   const std::vector<T> &new_row) const {
        assert(0 <= row_index && row_index < row_size);
        assert(static_cast<int>(new_row.size()) == column_size);
        std::vector<T> difference = new_row;
        const std::vector<T> current_row = get_row(row_index);
        for (int j = 0; j < column_size; ++j) {
            difference[j] -= current_row[j];
        }
        std::vector<T> column_vector(row_size, T());
        column_vector[row_index] = T(1);
        return rank_after_rank_one_update(column_vector, difference);
    }

    int rank_after_column_replacement(int column_index,
                                      const std::vector<T> &new_column) const {
        assert(0 <= column_index && column_index < column_size);
        assert(static_cast<int>(new_column.size()) == row_size);
        std::vector<T> difference = new_column;
        const std::vector<T> current_column = get_column(column_index);
        for (int i = 0; i < row_size; ++i) {
            difference[i] -= current_column[i];
        }
        std::vector<T> row_vector(column_size, T());
        row_vector[column_index] = T(1);
        return rank_after_rank_one_update(difference, row_vector);
    }

    int apply_rank_one_update(const std::vector<T> &column_vector,
                              const std::vector<T> &row_vector) {
        RankOneUpdateInfo info =
            analyze_rank_one_update(column_vector, row_vector);

        if (info.next_rank == matrix_rank + 1) {
            const int pivot_row = first_nonzero(info.column_residual);
            const int pivot_column = first_nonzero(info.row_residual);
            const std::vector<T> lambda = normalized_left_annihilator(
                pivot_row, info.column_residual[pivot_row]);
            const std::vector<T> rho = normalized_right_annihilator(
                pivot_column, info.row_residual[pivot_column]);

            const std::vector<std::vector<T>> old_row_space_right_inverse =
                row_space_right_inverse;

            append_column(column_space_basis, info.column_residual);
            column_space_left_inverse.push_back(lambda);

            row_space_basis.push_back(row_vector);
            for (int i = 0; i < matrix_rank; ++i) {
                for (int j = 0; j < column_size; ++j) {
                    row_space_basis[i][j] += info.alpha[i] * row_vector[j];
                }
            }

            std::vector<T> right_alpha(column_size, T());
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    right_alpha[i] +=
                        old_row_space_right_inverse[i][j] * info.alpha[j];
                }
            }

            row_space_right_inverse.assign(
                column_size, std::vector<T>(matrix_rank + 1, T()));
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    row_space_right_inverse[i][j] =
                        old_row_space_right_inverse[i][j] -
                        rho[i] * info.beta[j];
                }
                row_space_right_inverse[i][matrix_rank] =
                    T() - right_alpha[i] + rho[i] * info.schur;
            }
            ++matrix_rank;
            return matrix_rank;
        }

        if (!info.column_inside && info.row_inside) {
            const int pivot_row = first_nonzero(info.column_residual);
            const std::vector<T> lambda = normalized_left_annihilator(
                pivot_row, info.column_residual[pivot_row]);
            for (int i = 0; i < row_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    column_space_basis[i][j] += column_vector[i] * info.beta[j];
                }
            }
            std::vector<std::vector<T>> new_left_inverse =
                column_space_left_inverse;
            for (int i = 0; i < matrix_rank; ++i) {
                for (int j = 0; j < row_size; ++j) {
                    new_left_inverse[i][j] -= info.alpha[i] * lambda[j];
                }
            }
            column_space_left_inverse.swap(new_left_inverse);
            return matrix_rank;
        }

        if (info.column_inside && !info.row_inside) {
            const int pivot_column = first_nonzero(info.row_residual);
            const std::vector<T> rho = normalized_right_annihilator(
                pivot_column, info.row_residual[pivot_column]);
            for (int i = 0; i < matrix_rank; ++i) {
                for (int j = 0; j < column_size; ++j) {
                    row_space_basis[i][j] += info.alpha[i] * row_vector[j];
                }
            }
            std::vector<std::vector<T>> new_right_inverse =
                row_space_right_inverse;
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    new_right_inverse[i][j] -= rho[i] * info.beta[j];
                }
            }
            row_space_right_inverse.swap(new_right_inverse);
            return matrix_rank;
        }

        if (info.schur != T()) {
            std::vector<T> right_alpha(column_size, T());
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    right_alpha[i] +=
                        row_space_right_inverse[i][j] * info.alpha[j];
                }
            }
            for (int i = 0; i < matrix_rank; ++i) {
                for (int j = 0; j < column_size; ++j) {
                    row_space_basis[i][j] += info.alpha[i] * row_vector[j];
                }
            }
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    row_space_right_inverse[i][j] -=
                        right_alpha[i] * (info.beta[j] / info.schur);
                }
            }
            return matrix_rank;
        }

        const int removed = first_nonzero(info.alpha);
        const T removed_inverse = T(1) / info.alpha[removed];

        std::vector<T> right_alpha(column_size, T());
        for (int i = 0; i < column_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                right_alpha[i] += row_space_right_inverse[i][j] * info.alpha[j];
            }
        }

        std::vector<std::vector<T>> new_column_space_basis(
            row_size, std::vector<T>(matrix_rank - 1, T()));
        std::vector<std::vector<T>> new_row_space_basis(
            matrix_rank - 1, std::vector<T>(column_size, T()));
        std::vector<std::vector<T>> new_column_space_left_inverse(
            matrix_rank - 1, std::vector<T>(row_size, T()));
        std::vector<std::vector<T>> new_row_space_right_inverse(
            column_size, std::vector<T>(matrix_rank - 1, T()));

        int new_index = 0;
        for (int old = 0; old < matrix_rank; ++old) {
            if (old == removed) {
                continue;
            }
            for (int i = 0; i < row_size; ++i) {
                new_column_space_basis[i][new_index] =
                    column_space_basis[i][old] +
                    column_vector[i] * info.beta[old];
            }
            const T factor = info.alpha[old] * removed_inverse;
            for (int j = 0; j < column_size; ++j) {
                new_row_space_basis[new_index][j] =
                    row_space_basis[old][j] -
                    factor * row_space_basis[removed][j];
            }
            for (int j = 0; j < row_size; ++j) {
                new_column_space_left_inverse[new_index][j] =
                    column_space_left_inverse[old][j] -
                    factor * column_space_left_inverse[removed][j];
            }
            for (int i = 0; i < column_size; ++i) {
                new_row_space_right_inverse[i][new_index] =
                    row_space_right_inverse[i][old] +
                    right_alpha[i] * info.beta[old];
            }
            ++new_index;
        }

        column_space_basis.swap(new_column_space_basis);
        row_space_basis.swap(new_row_space_basis);
        column_space_left_inverse.swap(new_column_space_left_inverse);
        row_space_right_inverse.swap(new_row_space_right_inverse);
        --matrix_rank;
        return matrix_rank;
    }

    int apply_row_replacement(int row_index, const std::vector<T> &new_row) {
        assert(0 <= row_index && row_index < row_size);
        assert(static_cast<int>(new_row.size()) == column_size);
        std::vector<T> difference = new_row;
        const std::vector<T> current_row = get_row(row_index);
        for (int j = 0; j < column_size; ++j) {
            difference[j] -= current_row[j];
        }
        std::vector<T> column_vector(row_size, T());
        column_vector[row_index] = T(1);
        return apply_rank_one_update(column_vector, difference);
    }

    int apply_column_replacement(int column_index,
                                 const std::vector<T> &new_column) {
        assert(0 <= column_index && column_index < column_size);
        assert(static_cast<int>(new_column.size()) == row_size);
        std::vector<T> difference = new_column;
        const std::vector<T> current_column = get_column(column_index);
        for (int i = 0; i < row_size; ++i) {
            difference[i] -= current_column[i];
        }
        std::vector<T> row_vector(column_size, T());
        row_vector[column_index] = T(1);
        return apply_rank_one_update(difference, row_vector);
    }

  private:
    struct RankOneUpdateInfo {
        std::vector<T> alpha;
        std::vector<T> beta;
        std::vector<T> column_residual;
        std::vector<T> row_residual;
        bool column_inside = false;
        bool row_inside = false;
        T schur = T();
        int next_rank = 0;
    };

    RankOneUpdateInfo
    analyze_rank_one_update(const std::vector<T> &column_vector,
                            const std::vector<T> &row_vector) const {
        assert(static_cast<int>(column_vector.size()) == row_size);
        assert(static_cast<int>(row_vector.size()) == column_size);

        RankOneUpdateInfo info;
        info.alpha = multiply_left_inverse(column_vector);
        info.beta = multiply_right_inverse(row_vector);

        const std::vector<T> projected_column =
            reconstruct_column_vector(info.alpha);
        const std::vector<T> projected_row = reconstruct_row_vector(info.beta);

        info.column_residual.resize(row_size, T());
        info.row_residual.resize(column_size, T());
        info.column_inside = true;
        info.row_inside = true;
        for (int i = 0; i < row_size; ++i) {
            info.column_residual[i] = column_vector[i] - projected_column[i];
            if (info.column_residual[i] != T()) {
                info.column_inside = false;
            }
        }
        for (int j = 0; j < column_size; ++j) {
            info.row_residual[j] = row_vector[j] - projected_row[j];
            if (info.row_residual[j] != T()) {
                info.row_inside = false;
            }
        }

        info.schur = T(1);
        for (int i = 0; i < matrix_rank; ++i) {
            info.schur += info.beta[i] * info.alpha[i];
        }

        if (!info.column_inside && !info.row_inside) {
            info.next_rank = matrix_rank + 1;
        } else if (info.column_inside != info.row_inside) {
            info.next_rank = matrix_rank;
        } else {
            info.next_rank = info.schur == T() ? matrix_rank - 1 : matrix_rank;
        }
        return info;
    }

    std::vector<T>
    multiply_left_inverse(const std::vector<T> &column_vector) const {
        std::vector<T> result(matrix_rank, T());
        for (int i = 0; i < matrix_rank; ++i) {
            for (int j = 0; j < row_size; ++j) {
                result[i] += column_space_left_inverse[i][j] * column_vector[j];
            }
        }
        return result;
    }

    std::vector<T>
    multiply_right_inverse(const std::vector<T> &row_vector) const {
        std::vector<T> result(matrix_rank, T());
        for (int j = 0; j < column_size; ++j) {
            const T value = row_vector[j];
            if (value == T()) {
                continue;
            }
            for (int i = 0; i < matrix_rank; ++i) {
                result[i] += value * row_space_right_inverse[j][i];
            }
        }
        return result;
    }

    std::vector<T>
    reconstruct_column_vector(const std::vector<T> &coefficients) const {
        std::vector<T> result(row_size, T());
        for (int i = 0; i < row_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                result[i] += column_space_basis[i][j] * coefficients[j];
            }
        }
        return result;
    }

    std::vector<T>
    reconstruct_row_vector(const std::vector<T> &coefficients) const {
        std::vector<T> result(column_size, T());
        for (int i = 0; i < matrix_rank; ++i) {
            const T value = coefficients[i];
            if (value == T()) {
                continue;
            }
            for (int j = 0; j < column_size; ++j) {
                result[j] += value * row_space_basis[i][j];
            }
        }
        return result;
    }

    std::vector<T> normalized_left_annihilator(int pivot_row,
                                               const T &pivot_value) const {
        std::vector<T> result(row_size, T());
        result[pivot_row] = T(1);
        for (int j = 0; j < matrix_rank; ++j) {
            const T value = column_space_basis[pivot_row][j];
            if (value == T()) {
                continue;
            }
            for (int i = 0; i < row_size; ++i) {
                result[i] -= value * column_space_left_inverse[j][i];
            }
        }
        for (int i = 0; i < row_size; ++i) {
            result[i] /= pivot_value;
        }
        return result;
    }

    std::vector<T> normalized_right_annihilator(int pivot_column,
                                                const T &pivot_value) const {
        std::vector<T> result(column_size, T());
        result[pivot_column] = T(1);
        for (int i = 0; i < column_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                result[i] -= row_space_right_inverse[i][j] *
                             row_space_basis[j][pivot_column];
            }
        }
        for (int i = 0; i < column_size; ++i) {
            result[i] /= pivot_value;
        }
        return result;
    }

    static void append_column(std::vector<std::vector<T>> &matrix,
                              const std::vector<T> &column) {
        assert(static_cast<int>(matrix.size()) ==
               static_cast<int>(column.size()));
        for (int i = 0; i < static_cast<int>(matrix.size()); ++i) {
            matrix[i].push_back(column[i]);
        }
    }

    static int first_nonzero(const std::vector<T> &vector) {
        for (int i = 0; i < static_cast<int>(vector.size()); ++i) {
            if (vector[i] != T()) {
                return i;
            }
        }
        assert(false);
        return -1;
    }

    static std::vector<std::vector<T>>
    transpose_matrix(const std::vector<std::vector<T>> &matrix) {
        const int h = static_cast<int>(matrix.size());
        if (h == 0) {
            return {};
        }
        const int w = static_cast<int>(matrix[0].size());
        std::vector<std::vector<T>> transposed(w, std::vector<T>(h, T()));
        for (int i = 0; i < h; ++i) {
            assert(static_cast<int>(matrix[i].size()) == w);
            for (int j = 0; j < w; ++j) {
                transposed[j][i] = matrix[i][j];
            }
        }
        return transposed;
    }

    static std::vector<int>
    find_independent_columns(const std::vector<std::vector<T>> &matrix) {
        const int h = static_cast<int>(matrix.size());
        if (h == 0) {
            return {};
        }
        const int w = static_cast<int>(matrix[0].size());
        std::vector<std::vector<T>> b = matrix;
        int rank = 0;
        std::vector<int> pivot_columns;
        for (int column = 0; column < w; ++column) {
            int pivot = -1;
            for (int row = rank; row < h; ++row) {
                if (b[row][column] != T()) {
                    pivot = row;
                    break;
                }
            }
            if (pivot < 0) {
                continue;
            }
            if (pivot != rank) {
                b[pivot].swap(b[rank]);
            }

            const T inverse = T(1) / b[rank][column];
            for (int j = column; j < w; ++j) {
                b[rank][j] *= inverse;
            }
            for (int row = rank + 1; row < h; ++row) {
                const T factor = b[row][column];
                if (factor == T()) {
                    continue;
                }
                for (int j = column; j < w; ++j) {
                    b[row][j] -= b[rank][j] * factor;
                }
            }
            pivot_columns.push_back(column);
            ++rank;
            if (rank == h) {
                break;
            }
        }
        return pivot_columns;
    }

    static std::vector<std::vector<T>>
    inverse_matrix(std::vector<std::vector<T>> matrix) {
        const int n = static_cast<int>(matrix.size());
        for (int i = 0; i < n; ++i) {
            assert(static_cast<int>(matrix[i].size()) == n);
        }

        std::vector<std::vector<T>> inverse(n, std::vector<T>(n, T()));
        for (int i = 0; i < n; ++i) {
            inverse[i][i] = T(1);
        }

        for (int column = 0; column < n; ++column) {
            int pivot = -1;
            for (int row = column; row < n; ++row) {
                if (matrix[row][column] != T()) {
                    pivot = row;
                    break;
                }
            }
            assert(pivot >= 0);
            if (pivot != column) {
                matrix[pivot].swap(matrix[column]);
                inverse[pivot].swap(inverse[column]);
            }

            const T diagonal_inverse = T(1) / matrix[column][column];
            for (int j = 0; j < n; ++j) {
                matrix[column][j] *= diagonal_inverse;
                inverse[column][j] *= diagonal_inverse;
            }
            for (int row = 0; row < n; ++row) {
                if (row == column) {
                    continue;
                }
                const T factor = matrix[row][column];
                if (factor == T()) {
                    continue;
                }
                for (int j = 0; j < n; ++j) {
                    matrix[row][j] -= matrix[column][j] * factor;
                    inverse[row][j] -= inverse[column][j] * factor;
                }
            }
        }
        return inverse;
    }
};

#endif

#line 1 "math/matrix/dynamic-matrix-rank.hpp"



// 行列の階数を求め、外積 1 項更新後の階数を判定する。
// さらに、内部状態を O(k(r + c)) で更新しつつ変更後の階数を返せる。
// 現在の行列は左右の階数分解と片側逆元で保持する。
// 1 行差し替え・1 列差し替えはそれぞれ e_i (b-a_i)^T, (b-a_j) e_j^T に帰着する。
// 前処理は O(rc min(r, c) + k^2(r + c))、各更新・各判定は O(k(r + c)) である。
// T は体を成し、零判定と四則演算ができることを仮定する。

#include <cassert>
#include <vector>

template <class T> struct DynamicMatrixRank {
    int row_size = 0;
    int column_size = 0;
    int matrix_rank = 0;

    std::vector<std::vector<T>> column_space_basis;
    std::vector<std::vector<T>> row_space_basis;
    std::vector<std::vector<T>> column_space_left_inverse;
    std::vector<std::vector<T>> row_space_right_inverse;

    DynamicMatrixRank() = default;

    explicit DynamicMatrixRank(const std::vector<std::vector<T>> &matrix) {
        build(matrix);
    }

    void build(const std::vector<std::vector<T>> &matrix) {
        row_size = static_cast<int>(matrix.size());
        column_size = row_size == 0 ? 0 : static_cast<int>(matrix[0].size());
        for (int i = 1; i < row_size; ++i) {
            assert(static_cast<int>(matrix[i].size()) == column_size);
        }

        const std::vector<int> basis_columns = find_independent_columns(matrix);
        const std::vector<int> basis_rows =
            find_independent_columns(transpose_matrix(matrix));
        matrix_rank = static_cast<int>(basis_columns.size());
        assert(static_cast<int>(basis_rows.size()) == matrix_rank);

        std::vector<std::vector<T>> intersection_matrix(
            matrix_rank, std::vector<T>(matrix_rank, T()));
        for (int i = 0; i < matrix_rank; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                intersection_matrix[i][j] =
                    matrix[basis_rows[i]][basis_columns[j]];
            }
        }
        const std::vector<std::vector<T>> intersection_inverse =
            inverse_matrix(intersection_matrix);

        column_space_basis.assign(row_size, std::vector<T>(matrix_rank, T()));
        for (int i = 0; i < row_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                column_space_basis[i][j] = matrix[i][basis_columns[j]];
            }
        }

        row_space_basis.assign(matrix_rank, std::vector<T>(column_size, T()));
        for (int i = 0; i < matrix_rank; ++i) {
            for (int mid = 0; mid < matrix_rank; ++mid) {
                const T value = intersection_inverse[i][mid];
                if (value == T()) {
                    continue;
                }
                for (int j = 0; j < column_size; ++j) {
                    row_space_basis[i][j] += value * matrix[basis_rows[mid]][j];
                }
            }
        }

        column_space_left_inverse.assign(matrix_rank,
                                         std::vector<T>(row_size, T()));
        for (int i = 0; i < matrix_rank; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                column_space_left_inverse[i][basis_rows[j]] =
                    intersection_inverse[i][j];
            }
        }

        row_space_right_inverse.assign(column_size,
                                       std::vector<T>(matrix_rank, T()));
        for (int i = 0; i < matrix_rank; ++i) {
            row_space_right_inverse[basis_columns[i]][i] = T(1);
        }
    }

    void build() { build(materialize_matrix()); }

    int rank() const { return matrix_rank; }

    std::vector<T> get_row(int row_index) const {
        assert(0 <= row_index && row_index < row_size);
        std::vector<T> row(column_size, T());
        for (int i = 0; i < matrix_rank; ++i) {
            const T value = column_space_basis[row_index][i];
            if (value == T()) {
                continue;
            }
            for (int j = 0; j < column_size; ++j) {
                row[j] += value * row_space_basis[i][j];
            }
        }
        return row;
    }

    std::vector<T> get_column(int column_index) const {
        assert(0 <= column_index && column_index < column_size);
        std::vector<T> column(row_size, T());
        for (int i = 0; i < row_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                column[i] +=
                    column_space_basis[i][j] * row_space_basis[j][column_index];
            }
        }
        return column;
    }

    std::vector<std::vector<T>> materialize_matrix() const {
        std::vector<std::vector<T>> matrix(row_size,
                                           std::vector<T>(column_size, T()));
        for (int i = 0; i < row_size; ++i) {
            for (int mid = 0; mid < matrix_rank; ++mid) {
                const T value = column_space_basis[i][mid];
                if (value == T()) {
                    continue;
                }
                for (int j = 0; j < column_size; ++j) {
                    matrix[i][j] += value * row_space_basis[mid][j];
                }
            }
        }
        return matrix;
    }

    int rank_after_rank_one_update(const std::vector<T> &column_vector,
                                   const std::vector<T> &row_vector) const {
        return analyze_rank_one_update(column_vector, row_vector).next_rank;
    }

    int rank_after_row_replacement(int row_index,
                                   const std::vector<T> &new_row) const {
        assert(0 <= row_index && row_index < row_size);
        assert(static_cast<int>(new_row.size()) == column_size);
        std::vector<T> difference = new_row;
        const std::vector<T> current_row = get_row(row_index);
        for (int j = 0; j < column_size; ++j) {
            difference[j] -= current_row[j];
        }
        std::vector<T> column_vector(row_size, T());
        column_vector[row_index] = T(1);
        return rank_after_rank_one_update(column_vector, difference);
    }

    int rank_after_column_replacement(int column_index,
                                      const std::vector<T> &new_column) const {
        assert(0 <= column_index && column_index < column_size);
        assert(static_cast<int>(new_column.size()) == row_size);
        std::vector<T> difference = new_column;
        const std::vector<T> current_column = get_column(column_index);
        for (int i = 0; i < row_size; ++i) {
            difference[i] -= current_column[i];
        }
        std::vector<T> row_vector(column_size, T());
        row_vector[column_index] = T(1);
        return rank_after_rank_one_update(difference, row_vector);
    }

    int apply_rank_one_update(const std::vector<T> &column_vector,
                              const std::vector<T> &row_vector) {
        RankOneUpdateInfo info =
            analyze_rank_one_update(column_vector, row_vector);

        if (info.next_rank == matrix_rank + 1) {
            const int pivot_row = first_nonzero(info.column_residual);
            const int pivot_column = first_nonzero(info.row_residual);
            const std::vector<T> lambda = normalized_left_annihilator(
                pivot_row, info.column_residual[pivot_row]);
            const std::vector<T> rho = normalized_right_annihilator(
                pivot_column, info.row_residual[pivot_column]);

            const std::vector<std::vector<T>> old_row_space_right_inverse =
                row_space_right_inverse;

            append_column(column_space_basis, info.column_residual);
            column_space_left_inverse.push_back(lambda);

            row_space_basis.push_back(row_vector);
            for (int i = 0; i < matrix_rank; ++i) {
                for (int j = 0; j < column_size; ++j) {
                    row_space_basis[i][j] += info.alpha[i] * row_vector[j];
                }
            }

            std::vector<T> right_alpha(column_size, T());
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    right_alpha[i] +=
                        old_row_space_right_inverse[i][j] * info.alpha[j];
                }
            }

            row_space_right_inverse.assign(
                column_size, std::vector<T>(matrix_rank + 1, T()));
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    row_space_right_inverse[i][j] =
                        old_row_space_right_inverse[i][j] -
                        rho[i] * info.beta[j];
                }
                row_space_right_inverse[i][matrix_rank] =
                    T() - right_alpha[i] + rho[i] * info.schur;
            }
            ++matrix_rank;
            return matrix_rank;
        }

        if (!info.column_inside && info.row_inside) {
            const int pivot_row = first_nonzero(info.column_residual);
            const std::vector<T> lambda = normalized_left_annihilator(
                pivot_row, info.column_residual[pivot_row]);
            for (int i = 0; i < row_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    column_space_basis[i][j] += column_vector[i] * info.beta[j];
                }
            }
            std::vector<std::vector<T>> new_left_inverse =
                column_space_left_inverse;
            for (int i = 0; i < matrix_rank; ++i) {
                for (int j = 0; j < row_size; ++j) {
                    new_left_inverse[i][j] -= info.alpha[i] * lambda[j];
                }
            }
            column_space_left_inverse.swap(new_left_inverse);
            return matrix_rank;
        }

        if (info.column_inside && !info.row_inside) {
            const int pivot_column = first_nonzero(info.row_residual);
            const std::vector<T> rho = normalized_right_annihilator(
                pivot_column, info.row_residual[pivot_column]);
            for (int i = 0; i < matrix_rank; ++i) {
                for (int j = 0; j < column_size; ++j) {
                    row_space_basis[i][j] += info.alpha[i] * row_vector[j];
                }
            }
            std::vector<std::vector<T>> new_right_inverse =
                row_space_right_inverse;
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    new_right_inverse[i][j] -= rho[i] * info.beta[j];
                }
            }
            row_space_right_inverse.swap(new_right_inverse);
            return matrix_rank;
        }

        if (info.schur != T()) {
            std::vector<T> right_alpha(column_size, T());
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    right_alpha[i] +=
                        row_space_right_inverse[i][j] * info.alpha[j];
                }
            }
            for (int i = 0; i < matrix_rank; ++i) {
                for (int j = 0; j < column_size; ++j) {
                    row_space_basis[i][j] += info.alpha[i] * row_vector[j];
                }
            }
            for (int i = 0; i < column_size; ++i) {
                for (int j = 0; j < matrix_rank; ++j) {
                    row_space_right_inverse[i][j] -=
                        right_alpha[i] * (info.beta[j] / info.schur);
                }
            }
            return matrix_rank;
        }

        const int removed = first_nonzero(info.alpha);
        const T removed_inverse = T(1) / info.alpha[removed];

        std::vector<T> right_alpha(column_size, T());
        for (int i = 0; i < column_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                right_alpha[i] += row_space_right_inverse[i][j] * info.alpha[j];
            }
        }

        std::vector<std::vector<T>> new_column_space_basis(
            row_size, std::vector<T>(matrix_rank - 1, T()));
        std::vector<std::vector<T>> new_row_space_basis(
            matrix_rank - 1, std::vector<T>(column_size, T()));
        std::vector<std::vector<T>> new_column_space_left_inverse(
            matrix_rank - 1, std::vector<T>(row_size, T()));
        std::vector<std::vector<T>> new_row_space_right_inverse(
            column_size, std::vector<T>(matrix_rank - 1, T()));

        int new_index = 0;
        for (int old = 0; old < matrix_rank; ++old) {
            if (old == removed) {
                continue;
            }
            for (int i = 0; i < row_size; ++i) {
                new_column_space_basis[i][new_index] =
                    column_space_basis[i][old] +
                    column_vector[i] * info.beta[old];
            }
            const T factor = info.alpha[old] * removed_inverse;
            for (int j = 0; j < column_size; ++j) {
                new_row_space_basis[new_index][j] =
                    row_space_basis[old][j] -
                    factor * row_space_basis[removed][j];
            }
            for (int j = 0; j < row_size; ++j) {
                new_column_space_left_inverse[new_index][j] =
                    column_space_left_inverse[old][j] -
                    factor * column_space_left_inverse[removed][j];
            }
            for (int i = 0; i < column_size; ++i) {
                new_row_space_right_inverse[i][new_index] =
                    row_space_right_inverse[i][old] +
                    right_alpha[i] * info.beta[old];
            }
            ++new_index;
        }

        column_space_basis.swap(new_column_space_basis);
        row_space_basis.swap(new_row_space_basis);
        column_space_left_inverse.swap(new_column_space_left_inverse);
        row_space_right_inverse.swap(new_row_space_right_inverse);
        --matrix_rank;
        return matrix_rank;
    }

    int apply_row_replacement(int row_index, const std::vector<T> &new_row) {
        assert(0 <= row_index && row_index < row_size);
        assert(static_cast<int>(new_row.size()) == column_size);
        std::vector<T> difference = new_row;
        const std::vector<T> current_row = get_row(row_index);
        for (int j = 0; j < column_size; ++j) {
            difference[j] -= current_row[j];
        }
        std::vector<T> column_vector(row_size, T());
        column_vector[row_index] = T(1);
        return apply_rank_one_update(column_vector, difference);
    }

    int apply_column_replacement(int column_index,
                                 const std::vector<T> &new_column) {
        assert(0 <= column_index && column_index < column_size);
        assert(static_cast<int>(new_column.size()) == row_size);
        std::vector<T> difference = new_column;
        const std::vector<T> current_column = get_column(column_index);
        for (int i = 0; i < row_size; ++i) {
            difference[i] -= current_column[i];
        }
        std::vector<T> row_vector(column_size, T());
        row_vector[column_index] = T(1);
        return apply_rank_one_update(difference, row_vector);
    }

  private:
    struct RankOneUpdateInfo {
        std::vector<T> alpha;
        std::vector<T> beta;
        std::vector<T> column_residual;
        std::vector<T> row_residual;
        bool column_inside = false;
        bool row_inside = false;
        T schur = T();
        int next_rank = 0;
    };

    RankOneUpdateInfo
    analyze_rank_one_update(const std::vector<T> &column_vector,
                            const std::vector<T> &row_vector) const {
        assert(static_cast<int>(column_vector.size()) == row_size);
        assert(static_cast<int>(row_vector.size()) == column_size);

        RankOneUpdateInfo info;
        info.alpha = multiply_left_inverse(column_vector);
        info.beta = multiply_right_inverse(row_vector);

        const std::vector<T> projected_column =
            reconstruct_column_vector(info.alpha);
        const std::vector<T> projected_row = reconstruct_row_vector(info.beta);

        info.column_residual.resize(row_size, T());
        info.row_residual.resize(column_size, T());
        info.column_inside = true;
        info.row_inside = true;
        for (int i = 0; i < row_size; ++i) {
            info.column_residual[i] = column_vector[i] - projected_column[i];
            if (info.column_residual[i] != T()) {
                info.column_inside = false;
            }
        }
        for (int j = 0; j < column_size; ++j) {
            info.row_residual[j] = row_vector[j] - projected_row[j];
            if (info.row_residual[j] != T()) {
                info.row_inside = false;
            }
        }

        info.schur = T(1);
        for (int i = 0; i < matrix_rank; ++i) {
            info.schur += info.beta[i] * info.alpha[i];
        }

        if (!info.column_inside && !info.row_inside) {
            info.next_rank = matrix_rank + 1;
        } else if (info.column_inside != info.row_inside) {
            info.next_rank = matrix_rank;
        } else {
            info.next_rank = info.schur == T() ? matrix_rank - 1 : matrix_rank;
        }
        return info;
    }

    std::vector<T>
    multiply_left_inverse(const std::vector<T> &column_vector) const {
        std::vector<T> result(matrix_rank, T());
        for (int i = 0; i < matrix_rank; ++i) {
            for (int j = 0; j < row_size; ++j) {
                result[i] += column_space_left_inverse[i][j] * column_vector[j];
            }
        }
        return result;
    }

    std::vector<T>
    multiply_right_inverse(const std::vector<T> &row_vector) const {
        std::vector<T> result(matrix_rank, T());
        for (int j = 0; j < column_size; ++j) {
            const T value = row_vector[j];
            if (value == T()) {
                continue;
            }
            for (int i = 0; i < matrix_rank; ++i) {
                result[i] += value * row_space_right_inverse[j][i];
            }
        }
        return result;
    }

    std::vector<T>
    reconstruct_column_vector(const std::vector<T> &coefficients) const {
        std::vector<T> result(row_size, T());
        for (int i = 0; i < row_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                result[i] += column_space_basis[i][j] * coefficients[j];
            }
        }
        return result;
    }

    std::vector<T>
    reconstruct_row_vector(const std::vector<T> &coefficients) const {
        std::vector<T> result(column_size, T());
        for (int i = 0; i < matrix_rank; ++i) {
            const T value = coefficients[i];
            if (value == T()) {
                continue;
            }
            for (int j = 0; j < column_size; ++j) {
                result[j] += value * row_space_basis[i][j];
            }
        }
        return result;
    }

    std::vector<T> normalized_left_annihilator(int pivot_row,
                                               const T &pivot_value) const {
        std::vector<T> result(row_size, T());
        result[pivot_row] = T(1);
        for (int j = 0; j < matrix_rank; ++j) {
            const T value = column_space_basis[pivot_row][j];
            if (value == T()) {
                continue;
            }
            for (int i = 0; i < row_size; ++i) {
                result[i] -= value * column_space_left_inverse[j][i];
            }
        }
        for (int i = 0; i < row_size; ++i) {
            result[i] /= pivot_value;
        }
        return result;
    }

    std::vector<T> normalized_right_annihilator(int pivot_column,
                                                const T &pivot_value) const {
        std::vector<T> result(column_size, T());
        result[pivot_column] = T(1);
        for (int i = 0; i < column_size; ++i) {
            for (int j = 0; j < matrix_rank; ++j) {
                result[i] -= row_space_right_inverse[i][j] *
                             row_space_basis[j][pivot_column];
            }
        }
        for (int i = 0; i < column_size; ++i) {
            result[i] /= pivot_value;
        }
        return result;
    }

    static void append_column(std::vector<std::vector<T>> &matrix,
                              const std::vector<T> &column) {
        assert(static_cast<int>(matrix.size()) ==
               static_cast<int>(column.size()));
        for (int i = 0; i < static_cast<int>(matrix.size()); ++i) {
            matrix[i].push_back(column[i]);
        }
    }

    static int first_nonzero(const std::vector<T> &vector) {
        for (int i = 0; i < static_cast<int>(vector.size()); ++i) {
            if (vector[i] != T()) {
                return i;
            }
        }
        assert(false);
        return -1;
    }

    static std::vector<std::vector<T>>
    transpose_matrix(const std::vector<std::vector<T>> &matrix) {
        const int h = static_cast<int>(matrix.size());
        if (h == 0) {
            return {};
        }
        const int w = static_cast<int>(matrix[0].size());
        std::vector<std::vector<T>> transposed(w, std::vector<T>(h, T()));
        for (int i = 0; i < h; ++i) {
            assert(static_cast<int>(matrix[i].size()) == w);
            for (int j = 0; j < w; ++j) {
                transposed[j][i] = matrix[i][j];
            }
        }
        return transposed;
    }

    static std::vector<int>
    find_independent_columns(const std::vector<std::vector<T>> &matrix) {
        const int h = static_cast<int>(matrix.size());
        if (h == 0) {
            return {};
        }
        const int w = static_cast<int>(matrix[0].size());
        std::vector<std::vector<T>> b = matrix;
        int rank = 0;
        std::vector<int> pivot_columns;
        for (int column = 0; column < w; ++column) {
            int pivot = -1;
            for (int row = rank; row < h; ++row) {
                if (b[row][column] != T()) {
                    pivot = row;
                    break;
                }
            }
            if (pivot < 0) {
                continue;
            }
            if (pivot != rank) {
                b[pivot].swap(b[rank]);
            }

            const T inverse = T(1) / b[rank][column];
            for (int j = column; j < w; ++j) {
                b[rank][j] *= inverse;
            }
            for (int row = rank + 1; row < h; ++row) {
                const T factor = b[row][column];
                if (factor == T()) {
                    continue;
                }
                for (int j = column; j < w; ++j) {
                    b[row][j] -= b[rank][j] * factor;
                }
            }
            pivot_columns.push_back(column);
            ++rank;
            if (rank == h) {
                break;
            }
        }
        return pivot_columns;
    }

    static std::vector<std::vector<T>>
    inverse_matrix(std::vector<std::vector<T>> matrix) {
        const int n = static_cast<int>(matrix.size());
        for (int i = 0; i < n; ++i) {
            assert(static_cast<int>(matrix[i].size()) == n);
        }

        std::vector<std::vector<T>> inverse(n, std::vector<T>(n, T()));
        for (int i = 0; i < n; ++i) {
            inverse[i][i] = T(1);
        }

        for (int column = 0; column < n; ++column) {
            int pivot = -1;
            for (int row = column; row < n; ++row) {
                if (matrix[row][column] != T()) {
                    pivot = row;
                    break;
                }
            }
            assert(pivot >= 0);
            if (pivot != column) {
                matrix[pivot].swap(matrix[column]);
                inverse[pivot].swap(inverse[column]);
            }

            const T diagonal_inverse = T(1) / matrix[column][column];
            for (int j = 0; j < n; ++j) {
                matrix[column][j] *= diagonal_inverse;
                inverse[column][j] *= diagonal_inverse;
            }
            for (int row = 0; row < n; ++row) {
                if (row == column) {
                    continue;
                }
                const T factor = matrix[row][column];
                if (factor == T()) {
                    continue;
                }
                for (int j = 0; j < n; ++j) {
                    matrix[row][j] -= matrix[column][j] * factor;
                    inverse[row][j] -= inverse[column][j] * factor;
                }
            }
        }
        return inverse;
    }
};
Back to top page