18 #include "absl/strings/str_format.h"
31 const RowIndex num_rows = basis_matrix.
num_rows();
32 const ColIndex num_cols = basis_matrix.
num_cols();
38 basis_matrix_ = &basis_matrix;
41 lower_.
Reset(num_rows, num_cols);
42 upper_.
Reset(num_rows, num_cols);
43 permuted_lower_.
Reset(num_cols);
44 permuted_upper_.
Reset(num_cols);
45 permuted_lower_column_needs_solve_.
assign(num_cols,
false);
46 contains_only_singleton_columns_ =
true;
52 ExtractSingletonColumns(basis_matrix, row_perm, col_perm, &
index);
53 ExtractResidualSingletonColumns(basis_matrix, row_perm, col_perm, &
index);
54 int stats_num_pivots_without_fill_in =
index;
55 int stats_degree_two_pivot_columns = 0;
60 basis_matrix, *row_perm, *col_perm, &singleton_column_, &singleton_row_);
63 const int end_index =
std::min(num_rows.value(), num_cols.value());
65 parameters_.markowitz_singularity_threshold();
66 while (
index < end_index) {
75 const int64 min_markowitz = FindPivot(*row_perm, *col_perm, &pivot_row,
76 &pivot_col, &pivot_coefficient);
83 std::abs(pivot_coefficient) <= singularity_threshold) {
84 const std::string error_message = absl::StrFormat(
85 "The matrix is singular! pivot = %E", pivot_coefficient);
86 VLOG(1) <<
"ERROR_LU: " << error_message;
96 const int pivot_col_degree = residual_matrix_non_zero_.
ColDegree(pivot_col);
97 const int pivot_row_degree = residual_matrix_non_zero_.
RowDegree(pivot_row);
99 if (min_markowitz == 0) {
100 ++stats_num_pivots_without_fill_in;
101 if (pivot_col_degree == 1) {
102 RemoveRowFromResidualMatrix(pivot_row, pivot_col);
105 RemoveColumnFromResidualMatrix(pivot_row, pivot_col);
112 if (pivot_col_degree == 2) { ++stats_degree_two_pivot_columns; });
113 UpdateResidualMatrix(pivot_row, pivot_col);
116 if (contains_only_singleton_columns_) {
122 pivot_row, pivot_coefficient);
126 permuted_upper_.
column(pivot_col), pivot_row, pivot_coefficient);
131 (*col_perm)[pivot_col] = ColIndex(
index);
132 (*row_perm)[pivot_row] = RowIndex(
index);
136 stats_.pivots_without_fill_in_ratio.Add(
137 1.0 * stats_num_pivots_without_fill_in / end_index);
138 stats_.degree_two_pivot_columns.Add(1.0 * stats_degree_two_pivot_columns /
165 permuted_lower_.
Clear();
166 permuted_upper_.
Clear();
167 residual_matrix_non_zero_.
Clear();
168 col_by_degree_.
Clear();
169 examined_col_.clear();
170 is_col_by_degree_initialized_ =
false;
178 MatrixEntry(RowIndex r, ColIndex c,
Fractional coeff)
180 bool operator<(
const MatrixEntry& o)
const {
181 return (
row == o.row) ?
col < o.col :
row < o.row;
187 void Markowitz::ExtractSingletonColumns(
188 const CompactSparseMatrixView& basis_matrix,
RowPermutation* row_perm,
191 std::vector<MatrixEntry> singleton_entries;
192 const ColIndex num_cols = basis_matrix.num_cols();
193 for (ColIndex
col(0);
col < num_cols; ++
col) {
194 const ColumnView& column = basis_matrix.column(
col);
195 if (column.num_entries().value() == 1) {
196 singleton_entries.push_back(
197 MatrixEntry(column.GetFirstRow(),
col, column.GetFirstCoefficient()));
203 std::sort(singleton_entries.begin(), singleton_entries.end());
204 for (
const MatrixEntry e : singleton_entries) {
206 (*col_perm)[e.col] = ColIndex(*
index);
207 (*row_perm)[e.row] = RowIndex(*
index);
213 stats_.basis_singleton_column_ratio.Add(
static_cast<double>(*
index) /
217 bool Markowitz::IsResidualSingletonColumn(
const ColumnView& column,
220 int residual_degree = 0;
221 for (
const auto e : column) {
224 if (residual_degree > 1)
return false;
227 return residual_degree == 1;
230 void Markowitz::ExtractResidualSingletonColumns(
231 const CompactSparseMatrixView& basis_matrix,
RowPermutation* row_perm,
234 const ColIndex num_cols = basis_matrix.num_cols();
236 for (ColIndex
col(0);
col < num_cols; ++
col) {
238 const ColumnView& column = basis_matrix.column(
col);
239 if (!IsResidualSingletonColumn(column, *row_perm, &
row))
continue;
240 (*col_perm)[
col] = ColIndex(*
index);
241 (*row_perm)[
row] = RowIndex(*
index);
246 stats_.basis_residual_singleton_column_ratio.Add(
static_cast<double>(*
index) /
250 const SparseColumn& Markowitz::ComputeColumn(
const RowPermutation& row_perm,
265 if (permuted_lower_column_needs_solve_[
col]) {
269 const ColumnView&
input =
270 first_time ? basis_matrix_->
column(
col) : ColumnView(*lower_column);
273 permuted_lower_column_needs_solve_[
col] =
false;
274 return *lower_column;
280 if (lower_column->num_entries() == residual_matrix_non_zero_.
ColDegree(
col)) {
281 return *lower_column;
289 for (
const auto e : basis_matrix_->
column(
col)) {
290 lower_column->SetCoefficient(e.row(), e.coefficient());
293 lower_column->MoveTaggedEntriesTo(row_perm,
295 return *lower_column;
300 RowIndex* pivot_row, ColIndex* pivot_col,
305 while (!singleton_column_.empty()) {
306 const ColIndex
col = singleton_column_.back();
307 singleton_column_.pop_back();
316 if (residual_matrix_non_zero_.
ColDegree(
col) != 1)
continue;
321 if (contains_only_singleton_columns_) {
325 *pivot_row = e.row();
326 *pivot_coefficient = e.coefficient();
332 const SparseColumn& column = ComputeColumn(row_perm,
col);
333 if (column.IsEmpty())
continue;
335 *pivot_row = column.GetFirstRow();
336 *pivot_coefficient = column.GetFirstCoefficient();
339 contains_only_singleton_columns_ =
false;
344 while (!singleton_row_.empty()) {
345 const RowIndex
row = singleton_row_.back();
346 singleton_row_.pop_back();
354 if (residual_matrix_non_zero_.
RowDegree(
row) != 1)
continue;
358 const SparseColumn& column = ComputeColumn(row_perm,
col);
359 if (column.IsEmpty())
continue;
363 *pivot_coefficient = column.LookUpCoefficient(
row);
369 if (!is_col_by_degree_initialized_) {
370 is_col_by_degree_initialized_ =
true;
371 const ColIndex num_cols = col_perm.size();
372 col_by_degree_.
Reset(row_perm.size().value(), num_cols);
373 for (ColIndex
col(0);
col < num_cols; ++
col) {
375 const int degree = residual_matrix_non_zero_.
ColDegree(
col);
377 UpdateDegree(
col, degree);
384 examined_col_.clear();
385 const int num_columns_to_examine = parameters_.markowitz_zlatev_parameter();
386 const Fractional threshold = parameters_.lu_factorization_pivot_threshold();
387 while (examined_col_.size() < num_columns_to_examine) {
388 const ColIndex
col = col_by_degree_.
Pop();
391 const int col_degree = residual_matrix_non_zero_.
ColDegree(
col);
392 examined_col_.push_back(
col);
403 const int64 markowitz_lower_bound = col_degree - 1;
404 if (min_markowitz_number < markowitz_lower_bound)
break;
410 const SparseColumn& column = ComputeColumn(row_perm,
col);
411 DCHECK_EQ(column.num_entries(), col_degree);
415 max_magnitude =
std::max(max_magnitude, std::abs(e.coefficient()));
417 if (max_magnitude == 0.0) {
420 examined_col_.pop_back();
424 const Fractional skip_threshold = threshold * max_magnitude;
426 const Fractional magnitude = std::abs(e.coefficient());
427 if (magnitude < skip_threshold)
continue;
429 const int row_degree = residual_matrix_non_zero_.
RowDegree(e.row());
430 const int64 markowitz_number = (col_degree - 1) * (row_degree - 1);
432 if (markowitz_number < min_markowitz_number ||
433 ((markowitz_number == min_markowitz_number) &&
434 magnitude > std::abs(*pivot_coefficient))) {
435 min_markowitz_number = markowitz_number;
437 *pivot_row = e.row();
438 *pivot_coefficient = e.coefficient();
447 DCHECK_GE(min_markowitz_number, markowitz_lower_bound);
461 for (
const ColIndex
col : examined_col_) {
462 if (
col != *pivot_col) {
463 const int degree = residual_matrix_non_zero_.
ColDegree(
col);
467 return min_markowitz_number;
470 void Markowitz::UpdateDegree(ColIndex
col,
int degree) {
471 DCHECK(is_col_by_degree_initialized_);
483 singleton_column_.push_back(
col);
489 void Markowitz::RemoveRowFromResidualMatrix(RowIndex pivot_row,
490 ColIndex pivot_col) {
496 if (is_col_by_degree_initialized_) {
497 for (
const ColIndex
col : residual_matrix_non_zero_.
RowNonZero(pivot_row)) {
502 for (
const ColIndex
col : residual_matrix_non_zero_.
RowNonZero(pivot_row)) {
505 singleton_column_.push_back(
col);
511 void Markowitz::RemoveColumnFromResidualMatrix(RowIndex pivot_row,
512 ColIndex pivot_col) {
522 const RowIndex
row = e.row();
524 singleton_row_.push_back(
row);
529 void Markowitz::UpdateResidualMatrix(RowIndex pivot_row, ColIndex pivot_col) {
531 const SparseColumn& pivot_column = permuted_lower_.
column(pivot_col);
532 residual_matrix_non_zero_.
Update(pivot_row, pivot_col, pivot_column);
533 for (
const ColIndex
col : residual_matrix_non_zero_.
RowNonZero(pivot_row)) {
536 permuted_lower_column_needs_solve_[
col] =
true;
538 RemoveColumnFromResidualMatrix(pivot_row, pivot_col);
544 row_non_zero_.
clear();
545 deleted_columns_.
clear();
546 bool_scratchpad_.
clear();
547 num_non_deleted_columns_ = 0;
553 row_non_zero_.
clear();
554 row_non_zero_.
resize(num_rows.value());
555 deleted_columns_.
assign(num_cols,
false);
556 bool_scratchpad_.
assign(num_cols,
false);
557 num_non_deleted_columns_ = num_cols;
563 std::vector<RowIndex>* singleton_rows) {
564 const ColIndex num_cols = basis_matrix.
num_cols();
565 const RowIndex num_rows = basis_matrix.
num_rows();
568 Reset(num_rows, num_cols);
569 singleton_columns->clear();
570 singleton_rows->clear();
573 for (ColIndex
col(0);
col < num_cols; ++
col) {
575 deleted_columns_[
col] =
true;
576 --num_non_deleted_columns_;
580 ++row_degree_[e.row()];
585 for (RowIndex
row(0);
row < num_rows; ++
row) {
588 if (row_degree_[
row] == 1) singleton_rows->push_back(
row);
592 row_degree_[
row] = 0;
597 for (ColIndex
col(0);
col < num_cols; ++
col) {
599 int32 col_degree = 0;
601 const RowIndex
row = e.row();
607 col_degree_[
col] = col_degree;
608 if (col_degree == 1) singleton_columns->push_back(
col);
619 return --col_degree_[
col];
623 return --row_degree_[
row];
627 ColIndex pivot_col) {
628 DCHECK(!deleted_columns_[pivot_col]);
629 deleted_columns_[pivot_col] =
true;
630 --num_non_deleted_columns_;
633 row_degree_[pivot_row] = 0;
637 return deleted_columns_[
col];
641 auto& ref = row_non_zero_[
row];
643 const int end = ref.
size();
644 for (
int i = 0; i < end; ++i) {
645 const ColIndex
col = ref[i];
646 if (!deleted_columns_[
col]) {
647 ref[new_index] =
col;
651 ref.resize(new_index);
655 RowIndex
row)
const {
667 DCHECK(deleted_columns_[pivot_col]);
668 const int max_row_degree = num_non_deleted_columns_.value() + 1;
671 for (
const ColIndex
col : row_non_zero_[pivot_row]) {
673 bool_scratchpad_[
col] =
false;
681 const RowIndex
row = e.row();
682 if (
row == pivot_row)
continue;
687 if (e.coefficient() == 0.0 || row_degree_[
row] == max_row_degree)
continue;
696 const int kDeletionThreshold = 4;
697 if (row_non_zero_[
row].size() > row_degree_[
row] + kDeletionThreshold) {
702 MergeInto(pivot_row,
row);
710 MergeIntoSorted(pivot_row,
row);
715 void MatrixNonZeroPattern::MergeInto(RowIndex pivot_row, RowIndex
row) {
718 for (
const ColIndex
col : row_non_zero_[
row]) {
719 bool_scratchpad_[
col] =
true;
722 auto& non_zero = row_non_zero_[
row];
723 const int old_size = non_zero.
size();
724 for (
const ColIndex
col : row_non_zero_[pivot_row]) {
725 if (bool_scratchpad_[
col]) {
726 bool_scratchpad_[
col] =
false;
732 row_degree_[
row] += non_zero.
size() - old_size;
740 template <
typename V,
typename W>
741 void MergeSortedVectors(
const V& input_a, W* out) {
742 if (input_a.empty())
return;
743 const auto& input_b = *out;
744 int index_a = input_a.size() - 1;
745 int index_b = input_b.size() - 1;
746 int index_out = input_a.size() + input_b.size();
747 out->resize(index_out);
748 while (index_a >= 0) {
750 while (index_a >= 0) {
752 (*out)[index_out] = input_a[index_a];
758 if (input_a[index_a] > input_b[index_b]) {
759 (*out)[index_out] = input_a[index_a];
762 (*out)[index_out] = input_b[index_b];
773 void MatrixNonZeroPattern::MergeIntoSorted(RowIndex pivot_row, RowIndex
row) {
775 const auto&
input = row_non_zero_[pivot_row];
776 const auto& output = row_non_zero_[
row];
779 col_scratchpad_.resize(
input.size());
780 col_scratchpad_.resize(std::set_difference(
input.begin(),
input.end(),
781 output.begin(), output.end(),
782 col_scratchpad_.begin()) -
783 col_scratchpad_.begin());
786 for (
const ColIndex
col : col_scratchpad_) {
789 row_degree_[
row] += col_scratchpad_.
size();
790 MergeSortedVectors(col_scratchpad_, &row_non_zero_[
row]);
796 col_by_degree_.clear();
801 col_degree_.
assign(num_cols, 0);
802 col_index_.
assign(num_cols, -1);
803 col_by_degree_.resize(max_degree + 1);
804 min_degree_ = num_cols.value();
809 DCHECK_LT(degree, col_by_degree_.size());
810 const int32 old_degree = col_degree_[
col];
811 if (degree != old_degree) {
812 const int32 old_index = col_index_[
col];
813 if (old_index != -1) {
814 col_by_degree_[old_degree][old_index] = col_by_degree_[old_degree].back();
815 col_index_[col_by_degree_[old_degree].back()] = old_index;
816 col_by_degree_[old_degree].
pop_back();
819 col_index_[
col] = col_by_degree_[degree].size();
820 col_degree_[
col] = degree;
821 col_by_degree_[degree].push_back(
col);
822 min_degree_ =
std::min(min_degree_, degree);
824 col_index_[
col] = -1;
825 col_degree_[
col] = 0;
831 while (col_by_degree_[min_degree_].empty()) {
833 if (min_degree_ == col_by_degree_.size())
return kInvalidCol;
835 const ColIndex
col = col_by_degree_[min_degree_].back();
836 col_by_degree_[min_degree_].pop_back();
837 col_index_[
col] = -1;
838 col_degree_[
col] = 0;
843 mapping_.
assign(num_cols.value(), -1);
844 free_columns_.clear();
849 ColIndex
col)
const {
850 if (mapping_[
col] == -1)
return empty_column_;
851 return columns_[mapping_[
col]];
856 if (mapping_[
col] != -1)
return &columns_[mapping_[
col]];
858 if (free_columns_.empty()) {
859 new_col_index = columns_.size();
862 new_col_index = free_columns_.back();
863 free_columns_.pop_back();
865 mapping_[
col] = new_col_index;
866 return &columns_[new_col_index];
871 free_columns_.push_back(mapping_[
col]);
872 columns_[mapping_[
col]].Clear();
878 free_columns_.clear();