LARSCH
(convex/monge/larsch.hpp)
- View this file on GitHub
- Last update: 2026-07-07 18:38:36+09:00
- Include:
#include "convex/monge/larsch.hpp"
Overview
LARSCH finds row minima of an implicit lower-triangular totally monotone matrix
online. For row i, only columns 0 through i are candidates.
Unlike ordinary SMAWK, rows are produced in increasing order. Matrix entries in later columns may therefore depend on answers already computed for earlier rows. This supports dynamic programs such as
\[dp[i] = \min_{0 \le j < i}\left(dp[j] + cost(j,i)\right)\]in linear time when the resulting triangular transition matrix is totally monotone.
The algorithm is due to Larmore and Schieber: On-line dynamic programming with applications to the prediction of RNA secondary structure.
Sequential Matrix Interface
template <class T>
class Larsch;
Construct it with the number of rows and a callback:
m1une::convex::Larsch<long long> optimizer(
row_count,
[&](int row, int column) {
return matrix_value(row, column);
});
Then call next_argmin() exactly once for each row, in order.
| Method | Meaning |
|---|---|
size() |
Total row count. |
processed_rows() |
Number of answers already returned. |
finished() |
Whether every row was processed. |
next_argmin() |
Returns the leftmost minimum column of the next row. |
When evaluating entry (row, column), column <= row. More importantly,
entries using column j + 1 are not accessed before row j has been answered.
This is the online guarantee that permits callbacks to depend on earlier DP
values.
Monge DAG Shortest Paths
The convenience function
template <class Cost>
auto larsch_shortest_path(int vertex_count, Cost cost);
solves the complete forward DAG recurrence:
\[dist[0] = 0,\qquad dist[i] = \min_{0 \le j < i}(dist[j] + cost(j,i)).\]It returns LarschShortestPathResult<T>:
| Member | Meaning |
|---|---|
distance[i] |
Minimum distance from vertex 0 to i. |
parent[i] |
Leftmost predecessor attaining that distance, or -1 for vertex 0. |
Requirement
The lower-triangular matrix queried by the algorithm must be totally monotone for minima. For the shortest-path helper, the matrix
\[A(i,j) = dist[j] + cost(j,i+1),\qquad 0 \le j \le i\]must have that property. A Monge edge-cost function is the standard sufficient condition.
The implementation does not verify total monotonicity.
Complexity
Both interfaces use $O(N)$ callback evaluations, $O(N)$ time, and $O(N)$ memory, assuming each callback evaluation takes $O(1)$.
Example
#include "convex/monge/larsch.hpp"
#include <vector>
int main() {
std::vector<long long> position = {0, 2, 5, 9, 14};
auto result = m1une::convex::larsch_shortest_path(
int(position.size()),
[&](int from, int to) {
long long distance = position[to] - position[from];
return distance * distance + 3;
});
}
Required by
Convex All
(convex/all.hpp)
Verified with
Code
#ifndef M1UNE_CONVEX_MONGE_LARSCH_HPP
#define M1UNE_CONVEX_MONGE_LARSCH_HPP 1
#include <cassert>
#include <functional>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
namespace m1une {
namespace convex {
template <class T>
class Larsch {
struct ReduceColumn;
struct ReduceRow {
int size;
std::function<T(int, int)> value;
int current_row = 0;
int boundary = 0;
std::unique_ptr<ReduceColumn> recursive;
explicit ReduceRow(int size_) : size(size_) {
if (size / 2 != 0) recursive = std::make_unique<ReduceColumn>(size / 2);
}
void set_value(std::function<T(int, int)> value_) {
value = std::move(value_);
if (recursive) {
recursive->set_value(
[&](int row, int column) { return value(row * 2 + 1, column); });
}
}
int next_argmin() {
int row = current_row++;
if (row % 2 == 0) {
int previous = boundary;
int next = row + 1 == size ? size - 1 : recursive->next_argmin();
boundary = next;
int best = previous;
for (int column = previous + 1; column <= next; column++) {
if (value(row, column) < value(row, best)) best = column;
}
return best;
}
return value(row, boundary) <= value(row, row) ? boundary : row;
}
};
struct ReduceColumn {
int size;
std::function<T(int, int)> value;
int current_row = 0;
std::vector<int> columns;
ReduceRow recursive;
explicit ReduceColumn(int size_) : size(size_), recursive(size_) {}
void set_value(std::function<T(int, int)> value_) {
value = std::move(value_);
recursive.set_value(
[&](int row, int column) { return value(row, columns[column]); });
}
int next_argmin() {
int row = current_row++;
int first = row == 0 ? 0 : row * 2 - 1;
int last = row * 2;
for (int column = first; column <= last; column++) {
while (int(columns.size()) != row &&
value(int(columns.size()) - 1, columns.back()) >
value(int(columns.size()) - 1, column)) {
columns.pop_back();
}
if (int(columns.size()) != size) columns.push_back(column);
}
return columns[recursive.next_argmin()];
}
};
int _size;
int _processed = 0;
std::unique_ptr<ReduceRow> _base;
public:
template <class Value>
explicit Larsch(int size, Value value)
: _size(size), _base(std::make_unique<ReduceRow>(size)) {
assert(size >= 0);
_base->set_value(std::function<T(int, int)>(std::move(value)));
}
int size() const {
return _size;
}
int processed_rows() const {
return _processed;
}
bool finished() const {
return _processed == _size;
}
int next_argmin() {
assert(!finished());
_processed++;
return _base->next_argmin();
}
};
template <class T>
struct LarschShortestPathResult {
std::vector<T> distance;
std::vector<int> parent;
};
template <class Cost>
auto larsch_shortest_path(int vertex_count, Cost cost)
-> LarschShortestPathResult<
std::decay_t<std::invoke_result_t<Cost, int, int>>> {
using T = std::decay_t<std::invoke_result_t<Cost, int, int>>;
assert(vertex_count >= 0);
LarschShortestPathResult<T> result;
result.distance.assign(vertex_count, T());
result.parent.assign(vertex_count, -1);
if (vertex_count <= 1) return result;
Larsch<T> optimizer(vertex_count - 1, [&](int row, int column) {
return result.distance[column] + cost(column, row + 1);
});
for (int vertex = 1; vertex < vertex_count; vertex++) {
int parent = optimizer.next_argmin();
result.parent[vertex] = parent;
result.distance[vertex] = result.distance[parent] + cost(parent, vertex);
}
return result;
}
} // namespace convex
} // namespace m1une
#endif // M1UNE_CONVEX_MONGE_LARSCH_HPP#line 1 "convex/monge/larsch.hpp"
#include <cassert>
#include <functional>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
namespace m1une {
namespace convex {
template <class T>
class Larsch {
struct ReduceColumn;
struct ReduceRow {
int size;
std::function<T(int, int)> value;
int current_row = 0;
int boundary = 0;
std::unique_ptr<ReduceColumn> recursive;
explicit ReduceRow(int size_) : size(size_) {
if (size / 2 != 0) recursive = std::make_unique<ReduceColumn>(size / 2);
}
void set_value(std::function<T(int, int)> value_) {
value = std::move(value_);
if (recursive) {
recursive->set_value(
[&](int row, int column) { return value(row * 2 + 1, column); });
}
}
int next_argmin() {
int row = current_row++;
if (row % 2 == 0) {
int previous = boundary;
int next = row + 1 == size ? size - 1 : recursive->next_argmin();
boundary = next;
int best = previous;
for (int column = previous + 1; column <= next; column++) {
if (value(row, column) < value(row, best)) best = column;
}
return best;
}
return value(row, boundary) <= value(row, row) ? boundary : row;
}
};
struct ReduceColumn {
int size;
std::function<T(int, int)> value;
int current_row = 0;
std::vector<int> columns;
ReduceRow recursive;
explicit ReduceColumn(int size_) : size(size_), recursive(size_) {}
void set_value(std::function<T(int, int)> value_) {
value = std::move(value_);
recursive.set_value(
[&](int row, int column) { return value(row, columns[column]); });
}
int next_argmin() {
int row = current_row++;
int first = row == 0 ? 0 : row * 2 - 1;
int last = row * 2;
for (int column = first; column <= last; column++) {
while (int(columns.size()) != row &&
value(int(columns.size()) - 1, columns.back()) >
value(int(columns.size()) - 1, column)) {
columns.pop_back();
}
if (int(columns.size()) != size) columns.push_back(column);
}
return columns[recursive.next_argmin()];
}
};
int _size;
int _processed = 0;
std::unique_ptr<ReduceRow> _base;
public:
template <class Value>
explicit Larsch(int size, Value value)
: _size(size), _base(std::make_unique<ReduceRow>(size)) {
assert(size >= 0);
_base->set_value(std::function<T(int, int)>(std::move(value)));
}
int size() const {
return _size;
}
int processed_rows() const {
return _processed;
}
bool finished() const {
return _processed == _size;
}
int next_argmin() {
assert(!finished());
_processed++;
return _base->next_argmin();
}
};
template <class T>
struct LarschShortestPathResult {
std::vector<T> distance;
std::vector<int> parent;
};
template <class Cost>
auto larsch_shortest_path(int vertex_count, Cost cost)
-> LarschShortestPathResult<
std::decay_t<std::invoke_result_t<Cost, int, int>>> {
using T = std::decay_t<std::invoke_result_t<Cost, int, int>>;
assert(vertex_count >= 0);
LarschShortestPathResult<T> result;
result.distance.assign(vertex_count, T());
result.parent.assign(vertex_count, -1);
if (vertex_count <= 1) return result;
Larsch<T> optimizer(vertex_count - 1, [&](int row, int column) {
return result.distance[column] + cost(column, row + 1);
});
for (int vertex = 1; vertex < vertex_count; vertex++) {
int parent = optimizer.next_argmin();
result.parent[vertex] = parent;
result.distance[vertex] = result.distance[parent] + cost(parent, vertex);
}
return result;
}
} // namespace convex
} // namespace m1une