Symbolic-Numeric Integration of Univariate Expressions based on Sparse Regression

The main symbolic-numeric integration algorithm is summarized in Algorithm 1. Let $x$ be the independent variable. For this paper, we assume that the input function to be integrated, $f:\mathbb{C}\rightarrow\mathbb{C}$, is a univariate expression of $x$. Additionally, we assume that $f$ is well defined in a closed subset of the complex plane with only isolated poles.

The candidates generation algorithm (Section 3.2) produces a list of generator expressions, $G_{0},G_{1},G_{2},\cdots$ (by applying Equations 12, 13 and 15), where $G_{i}=\sum_{k}a_{k}\theta_{k}(x)$, $a_{k}\in\mathbb{C}$, $\theta_{k}:\mathbb{C}\rightarrow\mathbb{C}$, and $\theta_{k}$ does not have a constant leading coefficient. Each $G_{i}$ is converted to a set of candidates, $T_{i}=\{\theta_{k}\}$. According to Eq. 15, $T_{i}\subset T_{i+1}$. We can ignore the constant leading coefficients because the final coefficients are calculated by the numerical part of the algorithm. Therefore, we can express $\theta_{k}$s as equivalence classes of expressions where $f(x)\sim g(x)$ iff $f/g$ is a non-zero constant. Using this convention, we may write $(\sin x)^{\prime}=\cos x$ and $(\cos x)^{\prime}=\sin x$.

For each $G_{i}$ and the corresponding $T_{i}$, Algorithm 1 calls Algorithm 2 to find the anti-derivative of the input. If there is an acceptable solution, it returns; otherwise, Algorithm 1 fetches the next generator, $G_{i+1}$, and tests the resulting candidate set, $T_{i+1}$. If no solution is found after trying $G_{L-1}$, it returns with a failure message. $L$ is a parameter specified by the user.

Algorithm 2 generates $n$ test points $x_{i}$ in complex plane, where $n$ is the number of the candidates, to create an $n$-by-$n$ matrix $\mathbf{A}$ and an $n$-element vector $\mathbf{b}$ filled, respectively, with the values of the derivatives of the candidate terms and the input function at the test points.

As discussed above, a potential complication at this stage is that the rows of $\mathbf{A}$ may be linearly-dependent. We remedy this problem by finding the linearly-dependent rows and removing them from $\mathbf{A}$ and $\mathbf{b}$ (Section 4.3).

Using full-rank $\mathbf{A}$ and $\mathbf{b}$, we find $\mathbf{q}$ such that $\mathbf{A}\mathbf{q}=\mathbf{b}$. As is discussed in Section 4.4, we cast the problem as an optimization task and use sparse regression to calculate $\mathbf{q}$. Finally, we put everything together and generate $y=\sum_{j}q_{j}\theta_{j}$, the candidate anti-derivative of the input.

Candidate generation is the core of the Symbolic-Numeric integration algorithm. First, we describe an algorithm for a subset of expressions amenable to simple treatment and then expand to the general algorithm.

One key observation is that the form of the anti-derivative of some functions is similar to their derivative forms. For example, $\int\cos x\,dx=\sin x=-(\cos x)^{\prime}$. The underlying reason is that these functions can be defined in term of the polynomials of the exponential function and its inverse, i.e., $f\in\mathbb{C}[e^{x},e^{-x}]$, $(e^{x})^{\prime}=e^{x}$, and $(e^{-x})^{\prime}=-e^{-x}$; therefore $\mathbb{C}[e^{x},e^{-x}]$ is closed under differentiation and integration. The functions with this property include the non-negative integer powers of $\sin x$, $\cos x$, $\sinh x$, $\cosh x$, and $e^{x}$. Let’s define

where $p_{k}\in\mathbb{Z}^{+}$. For these class of functions, we can obtain the candidate terms with repetitive differentiation. For example, let $f(x)=e^{x}\sin x$. We start with $T=\emptyset$ ($T$ is the set of candidates), differentiate $f$, and add the resulting terms to $T$ (ignoring the constant coefficients):

Next, we need to differentiate each term in $T$ that was not previously processed. The first term, $e^{x}\cos x$, gives back $f$ and the second term is equal to $f$. Therefore, $T=\{e^{x}\cos x,e^{x}\sin x\}$ is the final answer, consistent with $\int e^{x}\sin x\,dx=(e^{x}\sin x-e^{x}\cos x)/2$. The main idea is repetitive differentiation of the terms until we run out of new terms to process. Essentially, we are generating a closure of the input expression based on differentiation.

Every $f\in\mathbb{C}[e^{x},e^{-x}]$ has a finite number of candidate terms; therefore, the preceding algorithm always terminates. We prove this for functions defined according to Eq. 6. After differentiation, each resulting term is also in the form of Eq. 6 and $p_{1}$ remains the same in all of them. In addition, $p_{2}+p_{3}$ remain the same (one power of $\sin x$ converts to $\cos x$ or vice versa). Similarly, $p_{4}+p_{5}$ is a constant. Therefore, the sum of the powers, $p=\sum_{k}p_{k}$, is a constant. Since for each $p$ there are only a finite number of possible expressions, we run out of new terms in a finite number of steps and the algorithm terminates.

In general, most integrands are not in $\mathbb{C}[e^{x},e^{-x}]$ and candidates cannot all be found by just repetitive differentiation; however, repetitive differentiation is still the backbone of the general algorithm. The essence of candidate generation is integration by parts,

followed by repetitive differentiation. This process or something equivalent resides at the core of most symbolic integration packages. However, when applying this algorithm to complex expressions, multiple problems arise. At each step, there can be more than one possible split of the integrand into $u$ and $v^{\prime}$. It is unknown a priori which split will result in the correct solution. Therefore, we need to explore all possibilities. Doing this naively can generate an exponential number of candidates. In fact, there is no guarantee that the process even terminates. Realizing that many intermediate expressions are shared among different paths, we can manage the number of candidates by keeping track of the intermediate results, similar to how memoization reduces the complexity of recursive algorithms. A top-down recursive algorithm with memoization usually has an equivalent bottom-up dynamic programming dual. This is the algorithmic basis of the approach discussed below.

A more systematic way to manage the intermediate expressions generated by repetitive integration by parts is provided by the continuous homotopy operators, which automate integration by parts and are useful in integrating exact expressions composed of multiple functions of independent variables and finding conservation laws in the system of differential equations [13, 14, 15]. Here, our goal is to borrow some of the machinery of the homotopy operators methodology to enhance the generation of candidate expressions. We continue with a brief and incomplete discussion of the homotopy operators technique.

Let $\mathbf{u}=(u^{1},u^{2},\cdots,u^{N})$ be a vector of $N$ dependent variables of $x$. Let $f(\mathbf{u})$ be a function of $\mathbf{u}$ and its partial derivatives with respect to $x$. We define $u_{0x}^{i}=u^{i}$, $u_{x}^{i}=\partial{u^{i}}/\partial{x}$ and $u^{i}_{kx}=\partial^{k}{u^{i}}/\partial{x}^{k}$. The total derivative operator is defined as

where $M^{j}$ is the maximum order of the partial derivatives of $u^{j}$ present in $f$. The key to the homotopy operators technique is a method to invert $\mathcal{D}_{x}f$, namely $\int f\,dx=\mathcal{D}_{x}^{-1}f=\mathcal{H}_{\mathbf{u}}f$, where

is the one-dimensional homotopy operator, $\lambda$ is a dummy integration variable, and the notation $X[\lambda\mathbf{u}]$ means replacing $u^{i}$ with $\lambda u^{i}$, $u_{x}^{i}$ with $\lambda u_{x}^{i}$, and the same for other variables and partial derivatives in $X$. Moreover, $\mathcal{I}_{u^{j}}$ is defined as

For the proof of these equations, refer to [13, 15]. However, we can make intuitive sense of these equations. Eq 11 has two main roles. First, it replaces $u_{kx}$ (removed by the partial differentiation operator) with $u_{(k-1)x}$. This action is similar to integration by part, and, in fact, its proof is based on integration by parts. Second, it performs repetitive differentiation (by the sub-expression containing the powers of $\mathcal{D}_{x}$), reminiscent of the method discussed above for functions in $\mathbb{C}[e^{x},e^{-x}]$. Finally, the definite integral in Eq 10 on variable $\lambda$ calculates the coefficients of the candidates.

In this paper, we use sparse regression to find the coefficients; therefore, we do not need Eq 10. Our focus is on Eq 11, but we cannot use it in its current form because the expected integrands are generally not in the form of $f(\mathbf{u})$. The key is to rewrite the integrand, $f(x)$, is a usable form,

$g_{i}(x)$ is a function that can be easily integrated (see Section 3.3), and $n_{i}\in\mathbb{Z}^{+}$. For example, if $f(x)=\sin^{2}(x^{2}-1)$, then $f=u_{1}^{2}$ for $u_{1}=\sin(v_{1})$ and $v_{1}=x^{2}-1$.

To find candidates, we follow the same process as in Eq 11 but customized for expressions conformant to Eq 12. However, we cannot simply apply Eq 11 because, for a general non-exact expression, $M^{j}$ may not be known or even bounded. Instead, we unroll the outer summation in Eq 11 to generate an ordered list of candidate generators, $G_{0},G_{1},G_{2},\cdots$. We generate $G_{0}$ by integrating each $u_{i}$ in turn. We integrate the first factor of $f$, i.e. $u_{1}^{n_{1}}$, by multiplying $f$ by $v_{1}^{\prime}/v_{1}^{\prime}$ and split $f$ into $u_{1}v_{1}^{\prime}$ and $f/(u_{1}v_{1}^{\prime})$. Ignoring the constant factors, the second part can be written as $(v_{1}^{\prime})^{-1}\partial{f}/\partial{u_{1}}$. Considering that $u=g(v)$, we have $\int uv^{\prime}\,dx=\int g(v)v^{\prime}\,dx=\int g(v)\,dv$ (remember that $g$ is chosen to be easy to integrate). Therefore,

In $\left(1+\int g_{i}(v)\,dv\right)$, 1 represents the constant of integration. For the term on the right, 1 is not strictly necessary but improves performance by moving terms from $G_{i>0}$ to $G_{0}$. Next, we generate $G_{1},G_{2},\cdots$ by repetitive differentiation similar to Eq 11. In addition, we introduce the powers of $x$ into the results by integration-by-parts, assuming $v^{\prime}=1$ and $v=x$ in Eq 8,

Putting all together,

Let’s apply the algorithm to $f(x)=\cot^{4}x$. Decomposing $f$ based on Eq 12, we have $u=\cot v$, $v=x$, and $f=u^{4}$. Considering that $\int\cot v\,dv=\log(\sin v)$, we get $G_{0}=1+\log(\sin x)+\cot^{3}x+\log(\sin x)\cot^{3}x$, which does not solve the integral. Next, we move to $G_{1}$ by using Eq 15. This time, we have

Feeding $f$ and the candidate set derived from $G_{1}$ to Algorithm 2, we solve the integral,

Converting integrands to a form compatible with Eq. 12 and finding $\int g_{i}(v)\,dv$ in Eq. 13 are facilitated by the term-rewriting and rule definition functionality of JuliaSymbolics, which includes a layered architectural design [16] to decouple term rewriting from other features, and leverages the Julia’s metaprogramming features that have been inherited from other languages in the Lisp family.

On the lowest level, TermInterface.jl and Metatheory.jl packages [17, 10, 16] provide an efficient expression rewriting engine and rule definition eDSL. Rules are written as regular Julia expressions, thus parsed into S-expressions and then efficiently compiled to callable functions that perform pattern matching and support advanced conditional selection [18, 10, 16]. JuliaSymbolics also supports the definition of equational (bidirectional) rewrite rules that can be executed through an e-graph rewriting engine (adopted from [19]), performing equality saturation.

We use JuliaSymbolics’s rewriting capabilities to find the integrals of $g(x)$. Nearly fifty rules are sufficient to cover elementary functions (exponential, logarithmic, trigonometric, hyperbolic, and their inverses). The following code snippets show a few example rules.

  @rule integrate(^(~x, -1)) => log(~x)
  @rule integrate(^(~x, ~k)) => ^(~x, ~k+1)
  @rule integrate(^(sin(~x), ~k::is_neg)) =>
         integrate(^(csc(~x), -~k))
  @rule integrate(csc(~x)) =>
         log(sin(~x)^-1 - cos(~x)*sin(~x)^-1)

Processing rational expressions is one of the first steps in symbolic integration with a long and mathematically rich history (e.g., Hermite’s and Harowitz’ methods). In this paper, to stay with the general theme of our method, we use a symbolic-numeric (hybrid) algorithm to factor rational expressions based on finding the roots of the denominator (which is assumed to be a polynomial). Our method is broadly similar to the one proposed in [20]; however, to find the coefficients of the terms, we use sparse regression method (see below) instead of the method of residues.

In line 2 of Algorithm 2, $n$ random number $x_{i}$ (test points) are generated in, $\mathbb{D}_{d}$, an open disk of radius $d$ centered at the origin ($d$ is a parameter provided by the user). Therefore, $\mathbb{D}_{d}$ should contain the majority of the unique poles of $f(x)$. Note that a function like $\sec(x)$ has an infinity number of poles and $\mathbb{D}_{d}$ should contain at least one of them. For the sparse regression to work properly, test points need to be near the poles of $f(x)$. We use a variant of the Newton-Raphson method to facilitate probing the poles. The standard Newton-Raphson method to find zeros of $f(x)$ is the iteration $x^{j+1}=x^{j}-f(x^{j})/f^{\prime}(x^{j})$, where superscript depicts an order index and not a power. The zero-finding algorithm can be transformed to a pole-finding one by changing the minus sign to plus. Therefore,

is used by Algorithm 2 on line 4 to move the test points toward the poles.

Matrix $\mathbf{A}$ in Algorithm 2 can be rank-deficient. In order to find acceptable regression solutions, we need to remove linearly-dependent columns from $\mathbf{A}$. This task is achieved by using pivoted QR algorithm [21]. In short, $\mathbf{A}$ is decomposed to $\mathbf{A}=\mathbf{Q}\mathbf{R}\mathbf{P}^{T}$, where $\mathbf{Q}$ is an orthogonal rotation matrix, $\mathbf{R}$ is an upper triangular matrix, and $\mathbf{P}$ is the permutation matrix. We find the columns linearly dependent on other columns by locating small values on the diagonal of $\mathbf{R}$ (here, small means having an absolute value less than a parameter $\epsilon$, usually taken to be $10^{-6}$) and convert $\mathbf{P}$ to the set $R$, used by Algorithm 2, that lists indices of linearly-dependent columns.

In Algorithm 2, we need to find $\mathbf{q}$ such that $\mathbf{A}\mathbf{q}=\mathbf{b}$. If $\mathbf{A}$ is low dimensional, this can be done simply as $\mathbf{q}=\mathbf{A}^{-1}\mathbf{b}$ (by construction, $\mathbf{A}$ is a square matrix). However, this process has the drawback of tending to use all the terms, even those with numerically small coefficients, which obscures the results and differs from the expected answer to integration problems. We need a parsimonious (sparse) model such that $\mathbf{q}$ has the minimum number of non-zero elements while still solves $\mathbf{A}\mathbf{q}=\mathbf{b}$ with an acceptable accuracy. We can achieve this by recasting the problem as an optimization problem to solve

for $\mathbf{q}$, where $i\in\{1,2~{}\}$, and $\lambda$ is a regularization parameter. In this paper, we use the sequential thresholded least-squares (STLSQ) algorithm, which is a component of Sparse Identification of Nonlinear Dynamics (SINDy) technique [22] and uses $\ell_{2}$-norm. This method has been chosen due to robust behaviour within the scope of many problems related to SINDy. However, other sparse regression algorithms, e.g., the least absolute shrinkage and selection operation (LASSO) using $\ell_{1}$-norm [23], Elastic-Net or SR3 are also possible candidates. The sparse regression code is provided by DataDrivenDiffEq.jl [11].