# 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 C^{3} library

The Compressed Continuous Computation (C^{3}) 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

- BLAS
- LAPACK
- 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

using different univariate approximation utilities available in

#### 3.1.1 Approximation in an orthonormal basis

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

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

We will now go through the steps of using univariate approximation. We begin with including the 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`x`

- an array of input locations , where 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 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 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|