Categories

# SMAWK in C++

I recently implemented kmeans1d—discussed in a prior post—for efficiently performing globally optimal 1D k-means clustering. The implementation utilizes the SMAWK algorithm (Aggarwal et al., 1987), which calculates `argmin(i)` for each row i of an arbitrary `n × m` totally monotone matrix, in `O(m(1 + lg(n/m)))`.

I’ve factored out my SMAWK C++ code into the example below. In general, SMAWK works with an implicitly defined matrix, utilizing a function that returns a value corresponding to an arbitrary position in the matrix. An explicitly defined matrix is used in the example for the purpose of illustration.

The program prints the column indices corresponding to the minimum element of each row in a totally monotone matrix. The matrix is from monge.pdf—a course document that I found online.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #include #include #include #include #include using namespace std; typedef unsigned long ulong; /* * Internal implementation of the SMAWK algorithm. */ template void _smawk( const vector& rows, const vector& cols, const function& lookup, vector* result) { // Recursion base case if (rows.size() == 0) return; // ******************************** // * REDUCE // ******************************** vector _cols; // Stack of surviving columns for (ulong col : cols) { while (true) { if (_cols.size() == 0) break; ulong row = rows[_cols.size() – 1]; if (lookup(row, col) >= lookup(row, _cols.back())) break; _cols.pop_back(); } if (_cols.size() < rows.size()) _cols.push_back(col); } // Call recursively on odd-indexed rows vector odd_rows; for (ulong i = 1; i < rows.size(); i += 2) { odd_rows.push_back(rows[i]); } _smawk(odd_rows, _cols, lookup, result); unordered_map col_idx_lookup; for (ulong idx = 0; idx < _cols.size(); ++idx) { col_idx_lookup[_cols[idx]] = idx; } // ******************************** // * INTERPOLATE // ******************************** // Fill-in even-indexed rows ulong start = 0; for (ulong r = 0; r < rows.size(); r += 2) { ulong row = rows[r]; ulong stop = _cols.size() – 1; if (r < rows.size() – 1) stop = col_idx_lookup[(*result)[rows[r + 1]]]; ulong argmin = _cols[start]; T min = lookup(row, argmin); for (ulong c = start + 1; c <= stop; ++c) { T value = lookup(row, _cols[c]); if (c == start || value < min) { argmin = _cols[c]; min = value; } } (*result)[row] = argmin; start = stop; } } /* * Interface for the SMAWK algorithm, for finding the minimum value * in each row of an implicitly-defined totally monotone matrix. */ template vector smawk( const ulong num_rows, const ulong num_cols, const function& lookup) { vector result; result.resize(num_rows); vector rows(num_rows); iota(begin(rows), end(rows), 0); vector cols(num_cols); iota(begin(cols), end(cols), 0); _smawk(rows, cols, lookup, &result); return result; } #define NUM_ROWS 9 #define NUM_COLS 18 // SMAWK works on implicitly defined matrices, utilizing a function // that returns a value as a function of matrix indices. // An explicitly defined matrix is used here for the purpose of // illustration. // The matrix is from: // http://web.cs.unlv.edu/larmore/Courses/CSC477/monge.pdf. double matrix[NUM_ROWS][NUM_COLS] = { { 25, 21, 13,10,20,13,19,35,37,41,58,66,82,99,124,133,156,178}, { 42, 35, 26,20,29,21,25,37,36,39,56,64,76,91,116,125,146,164}, { 57, 48, 35,28,33,24,28,40,37,37,54,61,72,83,107,113,131,146}, { 78, 65, 51,42,44,35,38,48,42,42,55,61,70,80,100,106,120,135}, { 90, 76, 58,48,49,39,42,48,39,35,47,51,56,63, 80, 86, 97,110}, {103, 85, 67,56,55,44,44,49,39,33,41,44,49,56, 71, 75, 84, 96}, {123,105, 86,75,73,59,57,62,51,44,50,52,55,59, 72, 74, 80, 92}, {142,123,100,86,82,65,61,62,50,43,47,45,46,46, 58, 59, 65, 73}, {151,130,104,88,80,59,52,49,37,29,29,24,23,20, 28, 25, 31, 39}, }; int main(void) { auto lookup = [](ulong i,ulong j) { return matrix[i][j]; }; vector argmin = smawk(NUM_ROWS, NUM_COLS, lookup); for (ulong i = 0; i < argmin.size(); ++i) { cout << argmin[i] << endl; } return 0; }
view raw smawk.cpp hosted with ❤ by GitHub

References

Aggarwal, Alok, Maria M. Klawe, Shlomo Moran, Peter Shor, and Robert Wilber. “Geometric Applications of a Matrix-Searching Algorithm.” Algorithmica 2, no. 1 (November 1, 1987): 195–208.