Mathematical Preliminaries
First, let's do some mathematical manipulation.
The equation you start with is:
$$
z-A|z|=b
$$
We would like to reduce this to something we know how to solve, like a standard linear system:
$$
Ax=b
$$
To do so, we introduce a diagonal matrix $S_z$ which we call the "signature matrix of $z$". Now, let the diagonal entries $S_i$ of the matrix have the value 1 if $z_i\ge0$ and the value -1 otherwise. That is:
$$
S_{ij}=
\begin{cases}
\phantom{-}0 & i\ne j \\
\phantom{-}1 & z_i \ge 0 \\
-1 & \textrm{otherwise}
\end{cases}
$$
Given the above, we have that
$$
S_z z = |z|
$$
This allows us to make the following moves
$$
\begin{align}
z-A|z|&=b \\
z-AS_z z&=b \\
(I-AS_z)z&=b
\end{align}
$$
The above can be solved using a standard linear solver... provided we know the value of $S_z$.
Thinking About Algorithms
As it turns out, there isn't an efficient way to figure this out because the problem is NP-Hard (see, for instance, this source). To get around this there's a slew of cunning heuristic algorithms that usually obtain results in reasonable timeframes (but there cannot be any general purpose algorithm that performs well for all cases unless P=NP).
While these algorithms would be handy for larger matrices, we won't use them here because you've said "that usually n≤5". This implies a brute force solution in which we'll solve up to 32 linear systems.
You might find this worrisome because you've also said
But I need very fast algorithms, so I cannot afford to solve 2^4 or 2^5 linear systems.
However, I think your aversion to linear solving is based on an incorrect cost model. You may be imagining that solving linear systems is expensive. You also may be imagining that there is some trick that lets you avoid the computation. The problem is, this trick isn't cheap, so it must save you a lot of time to be worth using.
To the first point, solving linear systems is not particularly expensive. It's one of the most fundamental operations in numeric computing and huge resources have been put into accelerating it. The Berkeley View is that dense linear algebra is a fundamental computational motif: in many cases if we can convert our problem to it we can declare victory and rest well in the knowledge that we have an efficient solution.
Your case is even more promising. The SSE and AVX instruction sets give processors the ability to perform up to 32 double-precision or 64 single-precision floating-point operations per clock cycle (in the case of AVX512). The effect of this is that, for a small matrix, rather than operating on a single matrix element at a time a good algorithm might process an entire row at once.
To the second point, the tricks you might use to avoid solving linear systems are not necessarily cheap. This paper seems to have a nice method for solving these systems and states that
The algorithm solved 100% of the test problems to an accuracy of $10^{−8}$ by solving an average of 2.77 systems of linear equations and linear programs per AVE.
As discussed here the cost of solving a linear program is $O(n^\omega \log(n/\delta))$ where $\delta$ is the relative accuracy. That is, linear programming has a time complexity equivalent to matrix multiplication (though the computational cost could be higher). Solving a linear system has a similar cost. Therefore, the algorithm quoted above costs about the same as ~6 linear system solves, nominally about ~5x faster than the brute force solution. But, remember, the brute force solution is simple and you have a small dataset; this means it has a good chance of outperforming a more complex algorithm.
The quoted algorithm goes on to state a wall-time of ~0.0215s per problem.
Results
Here I've implemented a brute force approach based on the mathematics we developed previously.
We compile with
g++ -O3 -std=c++17 -I/usr/include/eigen3 38832-solving-absolute-value-systems.cpp -Wall -pedantic -march=native
Note the use of -O3 to enable all optimizations and -march=native to allow the compiler to use the full power of the chip's instruction set.
Running this gives the following timing results on an Intel(R) Xeon(R) E-2176M CPU @ 2.70GHz:
╔═══════════╦═══════════╦════════════════╦═════════════════╦═════════════════════╦════════════╦════════════╗
║ Data Type ║ Dimension ║ Good Solutions ║ Full Solutions ║ Truncated Solutions ║ Time/Solve ║ Time/Solve ║
║ ║ ║ ║ Close to Actual ║ Close to Actual ║ Full ║ Truncated ║
╠═══════════╬═══════════╬════════════════╬═════════════════╬═════════════════════╬════════════╬════════════╣
║ float ║ 4x4 ║ 10k/10k ║ 6810/10k ║ 4054/10k ║ 4.4μs ║ 2.1μs ║
╠═══════════╬═══════════╬════════════════╬═════════════════╬═════════════════════╬════════════╬════════════╣
║ double ║ 4x4 ║ 10k/10k ║ 6691/10k ║ 4071/10k ║ 4.4μs ║ 2.4μs ║
╠═══════════╬═══════════╬════════════════╬═════════════════╬═════════════════════╬════════════╬════════════╣
║ float ║ 5x5 ║ 10k/10k ║ 5768/10k ║ 2836/10k ║ 12.8μs ║ 6.5μs ║
╠═══════════╬═══════════╬════════════════╬═════════════════╬═════════════════════╬════════════╬════════════╣
║ double ║ 5x5 ║ 10k/10k ║ 5710/10k ║ 2839/10k ║ 14.8μs ║ 7.6μs ║
╚═══════════╩═══════════╩════════════════╩═════════════════╩═════════════════════╩════════════╩════════════╝
The data type and size should be obvious. Good Solutions indicates that in all cases the program found a $z$ such that $z - A|z|$ was within $10^{-4}$ relative error of $b$. Close to Actual indicates how many of the returned $z$ values were within $10^{-4}$ of the actual $z$ value used to generate the problem. Full values are for instances in which the program was run for all values of $S_z$ while Truncated values are for instances in which the program made early exits when it found a solution that was good enough.
Early exits halve the solution time and produce results with low relative error, though these solutions may not correspond as often with the value used to generate the problem.
We don't always recover the $z$ used to generate the problem because of numeric issues and because multiple $z$ values can sometimes be good solutions.
Note that the worst-case timing result---14.8μs---is still 1453x faster than the timing results for the paper I quoted. This verifies our previous analysis: the simple algorithm out-performs the cunning algorithms, at least on datasets of this size. Even accounting for the smaller input sizes we've used and improved processor speeds since the quoted paper was written, the implementation we've developed here is still probably faster.
The Implementation
//Compile with: g++ -O3 -std=c++17 -I/usr/include/eigen3 38832-solving-absolute-value-systems.cpp -Wall -pedantic -march=native
#include <Eigen/Dense>
#include <bitset>
#include <chrono>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <random>
/// Solve the absolute value equation and return both the relative error and the
/// solution vector
template<class T, int DIM, bool early_exit>
std::pair<double, Eigen::Matrix<T, DIM, 1>> solve_absolute_value_equation(
const Eigen::Matrix<T, DIM, DIM>& A,
const Eigen::Matrix<T, DIM, 1>& b
){
// NOTE: Do not use `auto` with Eigen or many things will break
typedef Eigen::Matrix<T, DIM, DIM> m_type;
typedef Eigen::Matrix<T, DIM, 1> v_type;
double min_relative_error = std::numeric_limits<double>::infinity();
v_type best_z = v_type::Zero();
for(int signi=0; signi < (1<<DIM); signi++){ //Iterate over all Sz
// Construct Sz
const std::bitset<DIM> signs(signi);
Eigen::DiagonalMatrix<T, DIM> Sz;
for(int biti=0;biti<DIM;biti++){
Sz.diagonal()[biti] = signs.test(biti) ? 1.0 : -1.0;
}
// Solve the linear system
const m_type iasz = (m_type::Identity() - A * Sz);
const v_type z = iasz.colPivHouseholderQr().solve(b);
// Determine how close we are to the true solution
const v_type recovered = z - A * z.cwiseAbs();
const v_type diff = recovered - b;
const double relative_error = diff.norm() / b.norm(); // norm() is L2 norm
// Keep track of the best result
if(relative_error < min_relative_error){
min_relative_error = relative_error;
best_z = z;
}
// Return early if we find a promising solution (adjust value to your needs)
if(early_exit && relative_error < 1e-4){
break;
}
}
return std::make_pair(min_relative_error, best_z);
}
template<class T, int DIM, bool early_exit>
void timing_tests(){
// NOTE: Do not use `auto` with Eigen or many things will break
typedef Eigen::Matrix<T, DIM, DIM> m_type;
typedef Eigen::Matrix<T, DIM, 1> v_type;
// Get some timing statistics
constexpr auto NUM_TO_TEST = 10'000;
double total_time = 0;
int good_solutions = 0;
int close_to_actual = 0;
for(int i=0;i<NUM_TO_TEST;i++){
// Generate a problem instance with a known solution
const m_type A = m_type::Constant(0.5) - m_type::Random();
const v_type z_actual = v_type::Constant(0.5) - v_type::Random();
const v_type b = z_actual - A * z_actual.cwiseAbs();
const auto start = std::chrono::high_resolution_clock::now();
const auto [min_relative_error, best_z] = solve_absolute_value_equation<T, DIM, early_exit>(A, b);
const auto end = std::chrono::high_resolution_clock::now();
total_time += std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
const auto diff = best_z - z_actual;
const auto actual_rel_error = diff.norm() / z_actual.norm();
good_solutions += min_relative_error < 1e-4;
close_to_actual += actual_rel_error < 1e-4;
// std::cout << std::endl;
// std::cout << "Min relative error = " << min_relative_error << std::endl;
// std::cout << "Best z = " << best_z << std::endl;
// std::cout << "Actual z = " << z_actual << std::endl;
// std::cout << "Error versus the actual value = "<<actual_rel_error<<std::endl;
}
std::cout << "Average time per solution = " << (total_time / NUM_TO_TEST) << "us"<<std::endl;
std::cout << "Good solutions = " << good_solutions << " of " << NUM_TO_TEST << std::endl;
std::cout << "Solutions close to actual value = " << close_to_actual << " of " << NUM_TO_TEST << std::endl;
}
int main(){
// Set Eigen's random seed
srand(std::random_device()());
constexpr bool early_exit = false;
std::cout<<"float with dim=4"<<std::endl;
timing_tests<float, 4, early_exit>();
std::cout<<"double with dim=4"<<std::endl;
timing_tests<double, 4, early_exit>();
std::cout<<"float with dim=5"<<std::endl;
timing_tests<float, 5, early_exit>();
std::cout<<"double with dim=5"<<std::endl;
timing_tests<double, 5, early_exit>();
return 0;
}
```
