Compressed Continuous Computation

Table of Contents

1 Introduction

1.1 What is compressed?

When we refer to compression of an object we refer to the task of finding a representation of an object which takes up less space. Usually, we are able to reduce the amount of space representing an object by taking advantage of any structure present in the object. The purpose of this package is to compress functions into a form which both allows for efficient storage and efficient computation. We provide ways for converting a function (often represented using a function pointer) into the Function-Train format so that we may easily be able to perform computation with the functions in compressed format.

1.2 What is continuous computation?

Continuous computation refers to the philosophy of designing algorithms which work on continuous spaces instead of discrete arrays. Essentially, they relegate the "discretization" of functions to lower-levels of computation. Several benefits arise from such an approach:

  • Adaptation of discretizations can be performed more easily
  • Manipulation of functions is easier than manipulation of arrays of different sizes
  • Computation of continuous quantities such as derivatives or integrals can be obtained

To create this package we have been inspired by the work of Trefethen, Battles (2004), Townsend (2015) and others by the development of Chebfun, which performs computation with functions represented as Chebyshev polynomials.

1.3 C3 library

The Compressed Continuous Computation (C3) library is intended to make it easy to perform continuous linear and multilinear algebra with multidimensional functions. Common tasks include taking "matrix" decompositions of vector- or matrix-valued functions, adding or multiplying functions together, integrating multidimensional functions, and much much more.

Author: Alex A. Gorodetsky

Contact: alex@alexgorodetsky.com

Copyright: Massachusetts Institute of Technology 2014–2016, Alex A. Gorodetsky (2016)

Github: download License: BSD

2 Installation

The prerequisites for this package are

  1. BLAS
  2. LAPACK
  3. CMAKE

The installation uses CMAKE. I have tested it on Ubuntu and Mac OS X operating systems.

3 Tutorial

3.1 Univariate functions

We by approximating the function

\begin{equation*}
f(x) = x^2 + 3x + \cos(6x), \quad x \in [-2,3],
\end{equation*}

using different univariate approximation utilities available in $C^3$

3.1.1 Approximation in an orthonormal basis

The first form of the approximation that we consider is an expansion of Legendre polynomials:

\begin{equation*}
\hat{f}(x) = \sum_{i=1}^P a_i\phi_i, \quad \int_{-2}^3 \phi_i(x)\phi_j(x) dx = \delta_{i,j},
\end{equation*}

where we seek to find the optimal coefficients $a_i$ and the proper order $P$ to generate an accurate approximation.

We will now go through the steps of using univariate approximation. We begin with including the $C^3$ header and defining the function that we approximate

#include <stdlib.h>
#include "c3/c3.h"

int func(size_t N, const double * x, double * out, void * arg)
{
    double coeff = *(double *)arg;
    for (size_t ii = 0; ii < N; ii++){
        out[ii] = x[ii]*x[ii] + coeff*x[ii] + cos(6.0*x[ii]);
    }
    return 0;
}

The inputs of this function are

  • N - the number of locations at which to evaluate $f$
  • x - an array of input locations $N \times d$, where $d$ is the input dimension
  • out - an array at which to store the evaluations
  • arg - additional arguments, in our case the coefficient to the linear term

This function returns 0 if everything is successful and an integer otherwise.

Next, within the main() function we wrap this function into a $C^3$ specified struct:

// wrap the function
double coeff = 3.0;
struct Fwrap * fw = fwrap_create(1,"general-vec");
fwrap_set_fvec(fw,func,&coeff);

Next, we specify the approximation options and generate the approximation

// Create approximation options for orthogonal polynomial expansion (OPE)
// Legendre polynomials
struct OpeOpts * opts = ope_opts_alloc(LEGENDRE); 
// number of coefficients initially
ope_opts_set_start(opts,3); 
// stop adapting coefficients after they fall below tolerance
ope_opts_set_tol(opts,1e-15);
// number of coefficients to check
ope_opts_set_coeffs_check(opts,2); 
// lower bound of input (if not -1)
ope_opts_set_lb(opts,-2.0); 
// upper bound of input (if not 1)
ope_opts_set_ub(opts,3.0); 

// create approximation
struct GenericFunction * f = generic_function_approximate1d(POLYNOMIAL,opts,fw);

Thats it for the approximation portion! We now move on to computing with this function.

3.1.2 Approximation with piecewise polynomials

We now approximate the function with piecewise polynomials; each piece is approximated via an orthonormal basis.

The only difference between approximation with piecewise polynomials and an orthonormal basis is the choice of parameters. The final code shows the most commonly used parameters.

// Create approximation options for orthogonal polynomial expansion (OPE)
// Legendre polynomials, lower bound and upper bound
struct PwPolyOpts * opts = pw_poly_opts_alloc(LEGENDRE,-2.0,3.0);

pw_poly_opts_set_minsize(opts,1e-2); // set minimum size of pieces
pw_poly_opts_set_maxorder(opts,10); // maximum polynomial order of each piece
pw_poly_opts_set_tol(opts,1e-15); // threshold for coefficients of each piece
pw_poly_opts_set_coeffs_check(opts,2); // number of coefficients needed below tolerance

// create approximation
struct GenericFunction * f = generic_function_approximate1d(PIECEWISE,opts,fw);

Note the difference is that the first argument to generic_function_approximate1d is PIECEWISE instead of POLYNOMIAL.

3.1.3 Visualization of both approximations

3.1.4 Evaluation, differentiation, integration, and extremum seeking

We now assume that we have approximated a function and it is represented with the GenericFunction structure.

Evaluation of a univariate function requires two inputs: a GenericFunction struct and a double specifying the location of evaluation:

double pt   = 0.3;
double eval = generic_function_1d_eval(f,pt);

Differentiation of a univariate function requires two steps. First a new function is generated that represents the derivative. Next, the derivative function is evaluated.

double pt                   = 0.3;
struct GenericFunction * df = generic_function_deriv(f);
double deriv                = generic_function_1d_eval(df,pt);

Integration of a univariate function is simple:

double integral = generic_function_integral(f);

Finding the extremum of a function is a made by the following function call

double maxloc;
double absmax = generic_function_absmax(f,&maxvol,NULL);

The three inputs are:

  • f: the GenericFunction to maximize
  • maxvol: a reference to the location of the maximum that it finds
  • optargs: optional optimization arguments (in our case NULL)

The output is the absolute value of the extremum.

3.1.5 Complete example

Below is a complete working example for approximating a function in terms of a Legendre polynomial expansion. To compile it, add a link to the $C^3$ library with the -lc3 flag.

#include <stdlib.h>
#include "c3/c3.h"

int func(size_t N, const double * x, double * out, void * arg)
{
    double coeff = *(double *)arg;
    for (size_t ii = 0; ii < N; ii++){
        out[ii] = x[ii]*x[ii] + coeff*x[ii] + cos(6.0*x[ii]);
    }
    return 0;
}

int main()
{

  // wrap the function
  double coeff = 3.0;
  struct Fwrap * fw = fwrap_create(1,"general-vec");
  fwrap_set_fvec(fw,func,&coeff);


  // Create approximation options for legendre polynomial
  // orthogonal polynomial expansion options
  struct OpeOpts * opts = ope_opts_alloc(LEGENDRE); 
  // number of coefficients initially
  ope_opts_set_start(opts,3); 
  // stop adapting coefficients after they fall below tolerance
  ope_opts_set_tol(opts,1e-15);
  // number of coefficients to check
  ope_opts_set_coeffs_check(opts,2); 
  // lower bound of input (if not -1)
  ope_opts_set_lb(opts,-2.0); 
  // upper bound of input (if not 1)
  ope_opts_set_ub(opts,3.0); 

  // create approximation
  struct GenericFunction * f = generic_function_approximate1d(POLYNOMIAL,opts,fw);


  // Compare outputs
  // Evaluation
  double test_pt[1] = {0.3};
  double exact_eval[1],approx_eval[1];
  func(1,test_pt,exact_eval,&coeff);
  approx_eval[0] = generic_function_1d_eval(f,test_pt[0]);

  // Differentiation
  double exact_deriv = 2.0 * test_pt[0] + coeff - 6.0*sin(6.0*test_pt[0]);
  
  // Compute function that is the derivative
  struct GenericFunction * df = generic_function_deriv(f);
  double approx_deriv = generic_function_1d_eval(df,test_pt[0]);

  // Integration
  double exact_int = (pow(3.0,3)-pow(-2.0,3))/3.0 + 
                        coeff * (9.0-4.0)/2.0 + 
                        (sin(6.0*3)-sin(6.0*-2.0))/6.0;
  double approx_int = generic_function_integral(f);

  // Output
  printf("|Function   | Exact             | Approximate|\n");
  printf("|-----------+-------------------+------------------|\n");
  printf("|Evaluation | %3.15g | %3.15G|\n",exact_eval[0],approx_eval[0]);
  printf("|Derivative | %3.15g | %3.15G|\n",exact_deriv,approx_deriv);
  printf("|Integral   | %3.15g  | %3.15G|\n",exact_int,approx_int);


  // cleanup
  fwrap_destroy(fw);
  ope_opts_free(opts);
  generic_function_free(f);
  generic_function_free(df);
}
|Function   | Exact             | Approximate|
|-----------+-------------------+------------------|
|Evaluation | 0.762797905306913 | 0.762797905306967|
|Derivative | -2.24308578526917 | -2.24308578526956|
|Integral   | 18.9520733058713  | 18.9520733058714|

By . Share it—it's CC-BY-SA licensed.