Matrix and Vector L1 Norm

The L1 norm is defined for both vectors and matrices, we can easily write a C++ function to calculate it, but when possible it is better to use a more stable and generic implementation, as the one provided by the Boost Numeric uBLAS library.

L1 Norm for vectors

The L1 vector norm is defined as the sum of all absolute values of the vector elements. This is a simple L1 norm for double's vectors implementation:
namespace my
{
    double norm_1(const ublas::vector<double>& v)
    {
        auto adder = [](double partial, double value) { return partial + abs(value); }; // 1
        return std::accumulate(v.begin(), v.end(), 0.0, adder); // 2
    }
}
1. The natural implementation of this function is a one-liner, I have extract the lambda from the next line only for readability purpose. A simple adder is implemented, where the current partial result is passed in, increased with the absolute value of the current value, and returned.
2. The STL accumulate is used to sum all the element using our (1) custom adder. Note that the initial value should be explicitly set to double, since accumulate() deduces its return type from it.

L1 Norm for matrices

The L1 matrix norm is defined as the biggest L1 vector norm for each of its column. Here is a simple implementation:
namespace my
{
    double norm_1(const ublas::matrix<double>& m)
    {
        double norm = 0.0;
        for(unsigned int j = 0; j < m.size2(); ++j)
        {
            double current = 0.0; // 1
            for(unsigned int i = 0; i < m.size1(); ++i)
                current += abs(m(i,j));
            if(current > norm) // 2
                norm = current;
        }
        return norm;
    }
}
1. Calculate the L1 norm for a column in "current".
2. If the current candidate is bigger that the previously calculated candidate L1 norm, it is stored as new current L1 norm.

Examples

This piece of sample code shows how to calculate the L1 norm on a few vectors and matrices, using both my simple functions and the uBLAS versions:
// ...
namespace ublas = boost::numeric::ublas;

ublas::vector<double> v1(6);
ublas::vector<double> v2(6);
for(int i = 0; i < 6; ++i)
{
    v1[i] = i * (i%2 ? 1 : -1);
    v2[5 - i] = i;
}
ublas::vector<double> v3 = v1 * -1;

ublas::matrix<double> m1(3,3), m2(3,3);
for(int i = 0; i < 3; ++i)
{
    for(int j = 0; j < 3; ++j)
    {
        m1(i, j) = (i + j) * ((j+i)%2 ? 1 : -1);
        m2(2 - i, 2 - j) = i + j;
    }
}
ublas::matrix<double> m3 = m1 * -1;

std::cout << v1 << v2 << v3 << std::endl;
std::cout << m1 << std::endl << m2 << std::endl << m3 << std::endl;
std::cout << "Norm 1: " << ublas::norm_1(v1) << ' ' << my::norm_1(v1) << std::endl;
std::cout << "        " << ublas::norm_1(v2) << ' ' << my::norm_1(v2) << std::endl;
std::cout << "        " << ublas::norm_1(v3) << ' ' << my::norm_1(v3) << std::endl;
std::cout << "Norm 1: " << ublas::norm_1(m1) << ' ' << my::norm_1(m1) << std::endl;
std::cout << "        " << ublas::norm_1(m2) << ' ' << my::norm_1(m2) << std::endl;
std::cout << "        " << ublas::norm_1(m3) << ' ' << my::norm_1(m3) << std::endl;

No comments:

Post a Comment