Basic C++ for numerical computations: vectors

Norbert Pozar

This short note gives a few basic tips on how to efficiently write C++ code with a numerical performance in mind, focusing on vectors.

Introduction

Quick reminder before we begin: Your code should be as easy to understand as possible, easy to maintain and flexible. While it should be designed from the beginning with performance in mind, especially in terms of how data is stored in memory, flexibility is almost as important, especially in the early stages.

Vectors

C++ has the Standard Templates Library (STL) that provides some basic containers for storing elements. The vector STL container stores data in the same way the C array does, but with additional functionality.

After defining a variable of type vector, you can use it pretty much as a drop-in replacement of the C array. When optimization is enabled, the performance of vector and the C array is basically the same.

vector has a few advantages:

it remembers its size
it allocates and deallocates memory automatically, no need to touch pointers
fits nicely with other C++ libraries

Example: l_2-norm of a vector

Let us first look at the code for computing the l_2-norm of a vector using the standard C arrays.

#include <math>
int main() {
    int n = 100;
    double u[n];
    // initialize
    for (int i = 0; i < n; ++i) {
        u[i] = i;
    }
    // l2-norm
    double accum = 0.;
    for (int i = 0; i < n; ++i) {
        accum += u[i] * u[i];
    }
    double norm = sqrt(accum);
}

The same code can be easily rewritten using C++ and vector. Notice that we use u.size(), which returns the size of the vector. No need to remember the size anymore.

#include <vector>
#include <cmath>
using namespace std;
int main() {
    int n = 100;
    vector<double> u(n);  // create a vector with 100 elements, 0. by
                          // default
    // initialize (note that we use u.size() instead of n)
    for (int i = 0; i < u.size(); ++i) {
        u[i] = i;
    }
    // l2-norm
    double accum = 0.;
    for (int i = 0; i < n; ++i) {
        accum += u[i] * u[i];
    }
    double norm = sqrt(accum);
}

Of course, it makes sense to define a function that computes the l_2-norm. First the C code:

double l2_norm(double const* u, int n) {
    double accum = 0.;
    for (int i = 0; i < n; ++i) {
        accum += u[i] * u[i];
    }
    return sqrt(accum);
}

Because the C array does not know its size, we have to pass the size as the second argument of the function.

This becomes easier with vector:

double l2_norm(vector<double> const& u) {
    double accum = 0.;
    for (int i = 0; i < u.size(); ++i) {
        accum += u[i] * u[i];
    }
    return sqrt(accum);
}

The symbol & makes u a reference to a vector of doubles. This makes the code more efficient. Without &, u would be passed by value, which means that it would be copied into a local variable every time l2_norm is called. However, this is an unnecessary and slow operation. With &, if you call l2_norm in your code, for example:

vector<double> v;
...
double norm = l2_norm(v);

u inside l2_norm is just a different name for the same vector<double> that is stored in v. In particular, if you modify u inside l2_norm, v is modified as well. Hence the const keyword, which tells the caller that v will not be modified, even though it is passed by reference.

If your compiler supports the new C++11 standard (for g++ 4.6 and newer use the -std=c++0x command line option), the code can be simplified somewhat more using the new range-based loop:

double l2_norm(vector<double> const& u) {
    double accum = 0.;
    for (double x : u) {
        accum += x * x;
    }
    return sqrt(accum);
}

Example: Discretized Dirichlet energy

Let \Omega = (0,1). We want to find the minimizer of the Dirichlet energy

I(u):= \int_\Omega |u'(x)|^2 \;dx

over all functions in the set

\mathcal{K}:=\{v \in {\rm Lip}(\bar \Omega): v(0) = h_0, v(1) = h_1\}.

Here {\rm Lip} denotes all Lipschitz continuous functions on \bar\Omega. Lipschitz continuous functions are differentiable almost everywhere with bounded derivative, and therefore the integral is well-defined and finite.

To discretize this problem, we subdivide the domain into N \geq 1 equal subintervals (x_i, x_{i+1}), where x_i = i/N, i = 0, \ldots, N. Suppose that u is a piece-wise linear function on \Omega, linear on each (x_i, x_{i+1}). Function u is uniquely prescribed by its values v_i:= u(x_i) at the nodes x_i. We shall assume that the boundary values of u at x_0 = 0 and x_N = 1 are given values h_0, h_1. We introduce the discretized Dirichlet energy functional J(v_1, \ldots, v_{N-1}) as

J(v_1, \ldots, v_{N-1}):=\int_\Omega |u'(x)|^2 \;dx = \sum_{i=0}^{N-1}\frac{(v_{i+1}-v_i)^2}{\Delta x},

where \Delta x = x_{i+1} - x_i = 1/N.

Minimization of J over all vectors in \mathbb{R}^{N-1} is equivalent to minimization of the original energy I over all piece-wise linear functions on \Omega, linear on each subinterval (x_i, x_{i+1}), with boundary values u(0) = v_0 = h_0 and u(1) = v_N = h_1. It is a discretized version of the original problem in the sense that we restricted the search for the minimizer to a smaller set. Therefore the minimizer of the discretized problem will in general be only an approximation of the real minimizer.

An array or vector are a natural choice for storing the values of v_0, \ldots, v_N. The number of nodes (vector’s size) is N + 1.

It is easy to compute the discretized energy J, given the values v_0, \ldots, v_N.

double dirichlet_energy(vector<double> const& u) {
    int N = u.size() - 1;
    double accum = 0;
    for (int i = 0; i < N; ++i) {
        accum += pow(u[i+1] - u[i], 2.);
    }
    return N * accum;  // dx = 1 / N
}

The gradient of J has N-1 components:

\frac{\partial J}{\partial v_i} (v_1, \ldots, v_{N-1}) = -2 \frac{v_{i-1} - 2v_i + v_{i+1}}{\Delta x}.

However, it is convenient to store it as an array of N+1 elements and pad it on both sides by zeros.

vector<double> dirichlet_energy_gradient(vector<double> const& u) {
    int N = u.size() - 1;
    double c = -2. * N;  // c = -2. / dx
    vector<double> grad(u.size());  // vector of N+1 zeros
    for (int i = 1; i < N; ++i) {
        grad[i] = c * (u[i-1] - 2. * u[i] + u[i+1]);
    }
    return grad;
}

In a high-performance code one often wants to avoid unnecessary memory allocations and deallocations. Therefore it is better to return the gradient in a preallocated vector, which can be reused multiple times.

void dirichlet_energy_gradient(vector<double>& grad, vector<double> const& u) {
    // check that the vectors have the same size
    if (grad.size() != u.size()) {
        throw runtime_error("dirichlet_energy_gradient: size of vectors"
                "does not match");
    }
    int N = u.size() - 1;
    double c = -2. * N;  // c = -2. / dx
    for (int i = 1; i < N; ++i) {
        grad[i] = c * (u[i-1] - 2. * u[i] + u[i+1]);
    }
    grad[0] = grad[N] = 0.;  // do not forget the padding
}

We can also provide the function that returns the gradient directly, for convenience.

vector<double> dirichlet_energy_gradient(vector<double> const& u) {
    vector<double> grad(u.size());
    dirichlet_energy_gradient(grad, u);
    return grad;
}

Now we write a code to minimize the discretized Dirichlet energy using the steepest descent method (also called gradient descent), for simplicity. The minimization proceeds iteratively, at each iteration k we find the minimum of the functional J along the line of steepest descent, that is, along the line v_{k-1} - t g_{k-1}, where v_{k-1} is the solution after k-1 iterations, and g_{k-1} = \nabla J(v_{k-1}) is the gradient of J at v_{k-1}. The minimum t_k along this line must satisfy

\nabla J(v_{k-1} - t_k g_{k-1}) \cdot g_{k-1} = 0.

Since \nabla J is linear in this case, we can rewrite it as

0 = \nabla J(v_{k-1}) \cdot g_{k-1} - t_k \nabla J(g_{k-1}) \cdot g_{k-1} = g_{k-1} \cdot g_{k-1} - t_k \nabla J(g_{k-1}) \cdot g_{k-1}.

This gives a simple formula for t_k:

t_k = \frac{\|g_{k-1}\|^2}{\nabla J(g_{k-1}) \cdot g_{k-1}}.

The following full algorithm is implemented in find_dirichlet_minimizer below:

Set v_0 to be the initial guess. For k = 1, \ldots, k_{\rm max} do:

Compute the gradient, g_{k-1} = \nabla J(v_{k-1}).

If \|g_{k-1}\| \leq tol, exit. v_{k-1} is the resulting approximation and k-1 is the number of iterations performed.

Compute the minimizer along the line, t_k = \frac{\|g_{k-1}\|^2}{\nabla J(g_{k-1}) \cdot g_{k-1}}.

Set v_k = v_{k-1} - t_k g_{k-1}.

The code follows:

#include <cmath>
#include <vector>
#include <iostream>
#include <stdexcept>

using namespace std;

void check_size(vector<double> const& u, vector<double> const& v) {
    // check that the vectors have the same size
    if (u.size() != v.size()) {
        throw runtime_error("vector sizes do not match");
    }
}

double dirichlet_energy(vector<double> const& u) {
    int N = u.size() - 1;
    double accum = 0;
    for (int i = 0; i < N; ++i) {
        accum += pow(u[i+1] - u[i], 2.);
    }
    return N * accum;  // dx = 1 / N
}

void dirichlet_energy_gradient(vector<double>& grad, vector<double> const& u) {
    check_size(u, grad);
    int N = u.size() - 1;
    double c = -2. * N;  // c = -2. / dx
    for (int i = 1; i < N; ++i) {
        grad[i] = c * (u[i-1] - 2. * u[i] + u[i+1]);
    }
    grad[0] = grad[N] = 0.;  // do not forget the padding
}

double l2_norm(vector<double> const& u) {
    double accum = 0.;
    for (int i = 0; i < u.size(); ++i) {
        accum += u[i] * u[i];
    }
    return sqrt(accum);
}

double dot(vector<double> const& u, vector<double> const& v) {
    check_size(u, v);
    double accum = 0.;
    for (int i = 0; i < u.size(); ++i) {
        accum += u[i] * v[i];
    }
    return accum;
}

void add_to(vector<double>& u, double a, vector<double> const& v) {
    check_size(u, v);

    for (int i = 0; i < u.size(); ++i) {
        u[i] += a * v[i];
    }
}


int find_dirichlet_minimizer(vector<double>& u, double tolerance = 1e-5) {
    const int max_iter = 100000;
    vector<double> grad(u.size());
    vector<double> ggrad(u.size());

    int iter;
    for (iter = 0; iter < max_iter; ++iter) {
        dirichlet_energy_gradient(grad, u);
        double grad_norm = l2_norm(grad);

        if (grad_norm <= tolerance) {
            break;
        }

        // find the mimimum along the line u - t grad
        dirichlet_energy_gradient(ggrad, grad);
        double t = grad_norm * grad_norm / dot(ggrad, grad);

        // update
        add_to(u, -t, grad);
    }

    return iter;
}

int main() {
    // test the minimization
    int n = 128;

    vector<double> u(n + 1);

    for (int i = 0; i <= n; ++i) {
        double x = static_cast<double>(i) / n;
        u[i] = exp(x);
    }

    // quick test that the energy and gradient are compatible
    vector<double> grad(u.size());
    vector<double> v(u); // copy of u

    dirichlet_energy_gradient(grad, u);

    double delta = 1e-3;
    add_to(v, delta, grad);

    double Ju = dirichlet_energy(u);
    double Jv = dirichlet_energy(v);

    // should be small
    cout << "Relative Taylor difference = "
         << ((Ju + delta * dot(grad, grad) - Jv) / Ju) << endl;

    // find the minimizer
    int num_iter = find_dirichlet_minimizer(u);

    // report results
    cout << "Number of iterations = " << num_iter << endl;

    double J = dirichlet_energy(u);
    double J_exact = pow(u[n] - u[0], 2.);

    cout << "Energy at minimizer = " << J << endl;
    cout << "Difference from exact energy = " << (J - J_exact) << endl;

    dirichlet_energy_gradient(grad, u);
    cout << "Norm of the residual = " << l2_norm(grad) << endl;

    return 0;
}