expression compiler: automatic differentiation

[crowlogic]

Implementing Automatic Differentiator

This issue tracks the implementation of an automatic differentiator for the expression compiler. The compiler parses expressions into binary trees, and this differentiator will apply differentiation rules recursively based on the tree structure.

Algorithm Outline

1. Base Case:

   • For leaf nodes:
     • If a constant, derivative = 0.
     • If a variable, derivative = 1 for the differentiation variable, otherwise 0.

2. Differentiation Rules:

   • Addition/Subtraction Nodes: Derivative of a sum/difference is the sum/difference of the derivatives.
   • Multiplication Nodes (Product Rule): For a product $f * g$, derivative = $f' * g + f * g'$.
   • Division Nodes (Quotient Rule): For a division $f / g$, derivative = $(f' * g - f * g') / g^2$.
   • Exponential and Logarithmic Functions: Apply respective differentiation rules.
   • Trigonometric Functions: Apply differentiation rules for sine, cosine, etc.
   • Power Nodes: Apply power rule or generalized power rule.

3. Recursive Application:

   • Apply these rules to each child node to compute $f'$ and $g'$

4. Construct New Tree:

   • Construct a new tree representing the derivative.

5. Simplification (Optional):

   • Simplify the resulting tree to combine like terms, simplify constants, etc.

Stuff To Be Done And Whatnot

• Implement the base case logic for leaf nodes.
• Develop differentiation rules for each type of operation/function node.
• Create a function for recursive traversal and application of differentiation rules.
• Design and implement a tree constructor for the derivative expression.
• Add an optional simplification step post-differentiation.
• Thoroughly test the differentiator with various types of expressions.

[crowlogic]

package arb.functions.real;

import arb.Initializable;
import arb.Integer;
import arb.Real;
import arb.Typesettable;
import arb.documentation.BusinessSourceLicenseVersionOnePointOne;
import arb.documentation.TheArb4jLibrary;
import arb.expressions.nodes.DerivativeNode;
import junit.framework.TestCase;

/**
 * Decompiled {@link DerivativeNode} test function
 *
 * @see BusinessSourceLicenseVersionOnePointOne © terms of the
 *      {@link TheArb4jLibrary}
 */
public class TestCompiledDerivative implements
                                    RealFunctional<Object, RealFunction>,
                                    Typesettable,
                                    AutoCloseable,
                                    Initializable
{
  public boolean       isInitialized;
  public final Integer cℤ2     = new Integer("3");
  public final Integer cℤ1     = new Integer("2");
  public final Integer cℤ4     = new Integer("1");
  public final Integer cℤ3     = new Integer("0");
  public Real          a;
  public Real          b;
  public Real          c;
  public Real          ifuncℝ4 = new Real();
  public Real          ifuncℝ5 = new Real();
  public Integer       iℤ2     = new Integer();
  public Real          ifuncℝ6 = new Real();
  public Integer       iℤ1     = new Integer();
  public Real          ifuncℝ7 = new Real();
  public Real          ifuncℝ1 = new Real();
  public Real          ifuncℝ2 = new Real();
  public Real          ifuncℝ3 = new Real();
  public Real          ifuncℝ8 = new Real();

  public static void main(String args[])
  {
    try ( TestCompiledDerivative derivative = new TestCompiledDerivative())
    {
      derivative.a = Real.named("a").set(2);
      derivative.b = Real.named("b").set(4);
      derivative.c = Real.named("c").set(6);

      RealFunction d   = derivative.evaluate(null, 128);
      double       val = d.eval(2.3);
      TestCase.assertEquals(115.61999999999998, val);
      System.out.format("%s(2.3)=%s\n", d, val);
    }

  }

  @Override
  public Class<RealFunction> coDomainType()
  {
    return RealFunction.class;
  }

  @Override
  public RealFunction evaluate(Object in, int order, int bits, RealFunction result)
  {
    if (!isInitialized)
    {
      initialize();
    }
    RealFunction realFunction = new RealFunction()
    {

      @Override
      public Real evaluate(Real input, int order, int bits, Real res)
      {

        return a.add(b.mul(cℤ1.mul(input.pow(cℤ1.sub(cℤ4, bits, iℤ1), bits, ifuncℝ1), bits, ifuncℝ2), bits, ifuncℝ3),
                     bits,
                     ifuncℝ4)
                .add(c.mul(cℤ2.mul(input.pow(cℤ2.sub(cℤ4, bits, iℤ2), bits, ifuncℝ5), bits, ifuncℝ6), bits, ifuncℝ7),
                     bits,
                     ifuncℝ8);

      }

      @Override
      public String toString()
      {
        return TestCompiledDerivative.this.toString();
      }
    };
    return realFunction;
  }

  @Override
  public void initialize()
  {
    if (isInitialized)
    {
      throw new AssertionError("Already initialized");
    }
    else if (a == null)
    {
      throw new AssertionError("x-∂a*x+b*x²+c*x³⁄∂x.a is null");
    }
    else if (b == null)
    {
      throw new AssertionError("x-∂a*x+b*x²+c*x³⁄∂x.b is null");
    }
    else if (c == null)
    {
      throw new AssertionError("x-∂a*x+b*x²+c*x³⁄∂x.c is null");
    }
    else
    {
      isInitialized = true;
    }
  }

  @Override
  public void close()
  {
    cℤ2.close();
    cℤ1.close();
    cℤ4.close();
    cℤ3.close();
    ifuncℝ4.close();
    ifuncℝ5.close();
    iℤ2.close();
    ifuncℝ6.close();
    iℤ1.close();
    ifuncℝ7.close();
    ifuncℝ1.close();
    ifuncℝ2.close();
    ifuncℝ3.close();
    ifuncℝ8.close();
  }

  @Override
  public String toString()
  {
    return "x➔∂a*x+b*x²+c*x³/∂x";
  }

  @Override
  public String typeset()
  {
    return "a + b \\cdot 2 \\cdot {x}^{(\\left(2-1\\right))} + c \\cdot 3 \\cdot {x}^{(\\left(3-1\\right))}";
  }
}

[crowlogic]

To implement the differentiate method for the FunctionNode class, we need to apply the chain rule, as mentioned earlier. Here's how you can implement this method based on the information provided in the FunctionNode class. This implementation considers both built-in and contextual functions, as well as the requirement that a function may be specialized or generically invoked:

Key Points for Implementation

1. Differentiate the Argument: Apply the differentiation recursively on the argument node of the function.

2. Differentiate the Function Itself: Determine the derivative of the function with respect to its argument. This may involve either looking up a predefined derivative or implementing the derivative calculation directly.

3. Multiply the Results: Combine the results using the chain rule: $$ f'(g(x)) \cdot g'(x) $$.

Implementing the Differentiate Method

Here's an example implementation to guide you:

@Override
public Node<D, R, F> differentiate(VariableNode<D, R, F> variable) {
    // Step 1: Differentiate the argument (g'(x)).
    Node<D, R, F> argDerivative = arg.differentiate(variable);

    // Step 2: Differentiate the function (f'(g(x))).
    Node<D, R, F> functionDerivative = differentiateFunction();

    // Step 3: Apply the chain rule: f'(g(x)) * g'(x).
    return new MultiplicationNode<>(expression, functionDerivative, argDerivative);
}

/**
 * Returns the node representing the derivative of the function.
 * This will vary based on whether the function is built-in or contextual.
 */
private Node<D, R, F> differentiateFunction() {
    // Check if the function is built-in or contextual.
    if (isBuiltin()) {
        return differentiateBuiltinFunction();
    } else if (contextual) {
        return differentiateContextualFunction();
    } else {
        throw new UnsupportedOperationException("Cannot differentiate function: " + functionName);
    }
}

/**
 * Handles differentiation for built-in functions.
 */
private Node<D, R, F> differentiateBuiltinFunction() {
    switch (functionName) {
        case "sin":
            return new FunctionNode<>("cos", arg, expression); // derivative of sin is cos
        case "cos":
            return new NegationNode<>(expression, new FunctionNode<>("sin", arg, expression)); // derivative of cos is -sin
        case "exp":
            return this; // derivative of exp is exp
        // Add other built-in function derivatives
        default:
            throw new UnsupportedOperationException("Derivative not implemented for function: " + functionName);
    }
}

/**
 * Handles differentiation for contextual functions.
 */
private Node<D, R, F> differentiateContextualFunction() {
    // Add logic for differentiating contextual functions, potentially involving more advanced logic
    // depending on how contextual functions are defined and used in your framework.

    // As an example, you might retrieve a derivative function mapping if available:
    if (mapping != null && mapping.derivativeMapping != null) {
        return new FunctionNode<>(mapping.derivativeFunctionName, arg, expression);
    }
    // Alternatively, handle differentiation based on known properties of the function
    throw new UnsupportedOperationException("Contextual function differentiation not implemented: " + functionName);
}

Considerations

• Built-in Function Differentiation: For built-in functions, you can either hard-code the derivatives as shown or look them up if your framework supports such a table.

• Contextual Functions: For contextual functions, you may need to obtain or define their derivatives dynamically, depending on how they are registered and invoked within your framework.

• Error Handling: Handle functions that do not have derivatives or when the differentiation is not supported by throwing exceptions or providing meaningful feedback.

This implementation assumes that you have proper classes and methods like MultiplicationNode and NegationNode, and that your framework supports distinguishing between built-in and contextual functions. Adjust the logic as needed to fit the design and capabilities of your symbolic computation system.

Citations:
[1] https://ppl-ai-file-upload.s3.amazonaws.com/web/direct-files/37918610/8a0ff5f4-5a47-4a4f-b6f2-b871c663c989/paste.txt
[2] https://ppl-ai-file-upload.s3.amazonaws.com/web/direct-files/37918610/f4f9b1e9-fc44-4caa-8883-1fc0db30d2c5/paste-2.txt
[3] https://ppl-ai-file-upload.s3.amazonaws.com/web/direct-files/37918610/71dbc71c-3e90-4971-b110-f586f228daeb/paste-3.txt

[crowlogic]

For an n-ary product (product of n functions), the derivative follows this formula:

$$\frac{d}{dx}\left[\prod_{i=1}^{k}f_{i}(x)\right]=\sum_{i=1}^{k}\left(\left(\frac{d}{dx}f_{i}(x)\right)\prod_{j=1,j\neq i}^{k}f_{j}(x)\right) $$

Example for Three Functions

For three functions u, v, and w, the derivative is:

$$\frac{d(uvw)}{dx}=\frac{du}{dx}vw+u\frac{dv}{dx}w+uv\frac{dw}{dx} $$

Pattern

Each term in the sum is formed by:

1. Taking the derivative of one function
2. Multiplying it by all other functions unchanged
3. Adding all possible combinations

The rule extends to any number of functions following this same pattern.

Citations:
[1] https://www.mathcentre.ac.uk/resources/uploaded/mc-ty-product-2009-1.pdf
[2] https://en.wikipedia.org/wiki/Product_rule
[3] https://calcworkshop.com/derivatives/product-rule/
[4] https://www.cuemath.com/calculus/product-rule/
[5] https://byjus.com/maths/product-rule/

The derivative of a sum follows the linearity property of derivatives - you can differentiate each term separately and then sum the results:

$$\frac{d}{dx}\sum_{i=1}^{n}f_i(x) = \sum_{i=1}^{n}\frac{d}{dx}f_i(x) $$

Pattern

The derivative operator can move inside the summation because:

1. Differentiation is a linear operation
2. Each term is differentiated independently
3. The sum of derivatives equals the derivative of the sum

For example, if you have:

$$\frac{d}{dx}(f_1(x) + f_2(x) + ... + f_n(x)) = \frac{d}{dx}f_1(x) + \frac{d}{dx}f_2(x) + .. + \frac{d}{dx}f_n(x)$$

This is much simpler than the product rule because addition is a linear operation.

Citations:
[1] https://www.wolframalpha.com/input?input=d%2Fdx%28%E2%88%91_%7Bi%3D1%7D%5E%7Bn%7D+f%28i%29%29

[crowlogic]

For n-ary products, there is indeed a closed form for integration. Using the multinomial theorem and integration by parts, we can express it as:

$$ \int \prod_{i=1}^{n} f_i(x)dx = \sum_{k=1}^n \left((-1)^{k+1} \sum_{i_1 + ... + i_n = k} \int x^k \frac{f_1^{(i_1)}...f_n^{(i_n)}}{i_1!...i_n!} dx\right) $$

Where:

• $f_i^{(j)}$ represents the jth derivative of function i
• The inner sum is over all possible combinations where the indices sum to k
• The formula assumes all functions are infinitely differentiable

This is significantly more complex than differentiation, but it does provide a systematic way to compute the integral of a product of n functions.

Citations:
[1] https://www.youtube.com/watch?v=zObGec7JM3E
[2] https://en.wikipedia.org/wiki/Solution_in_closed_form

[crowlogic]

  public void testArcSinDerivative()
  {
    var f  = RealFunction.parse("∂arcsin(x)/∂x");
    var df = RealFunction.parse("1/√(1-x^2)");
    assertEquals(df.rootNode.toString(), f.rootNode.toString());
  }

[crowlogic]

Add support for f'(t) and f''(t) syntax for 1st and 2nd derivative, ...

Let's write it purely in terms of J_n(t).

The Bessel differential equation in standard form:

$$t^2J_n''(t) + tJ_n'(t) + (t^2 - n^2)J_n(t) = 0$$

To transform this into a Riccati equation, we use the substitution:

$$w(t) = \frac{J_n'(t)}{J_n(t)}$$

The transformation steps:

1. From the substitution:

$$J_n'(t) = w(t)J_n(t)$$

2. Differentiating:

$$J_n''(t) = w'(t)J_n(t) + w(t)^2J_n(t)$$

3. Substituting into the original equation:

$$t^2(w'(t)J_n(t) + w(t)^2J_n(t)) + t(w(t)J_n(t)) + (t^2 - n^2)J_n(t) = 0$$

4. Dividing by J_n(t):

$$t^2w'(t) + t^2w(t)^2 + tw(t) + (t^2 - n^2) = 0$$

5. The resulting Riccati equation:

$$w'(t) = -w(t)^2 - \frac{w(t)}{t} - (1 - \frac{n^2}{t^2})$$

This shows that w(t), the logarithmic derivative of the Bessel function J_n(t), satisfies this first-order nonlinear differential equation.

Ill be damned, it makes sense now

Yes, the logarithmic derivative w(t) = J_n'(t)/J_n(t) of the Bessel function J_n(t) satisfies the Riccati equation:

$$w'(t) = -w(t)^2 - \frac{w(t)}{t} - (1 - \frac{n^2}{t^2})$$

This comes directly from the fact that J_n(t) satisfies the Bessel differential equation:

$$t^2J_n''(t) + tJ_n'(t) + (t^2 - n^2)J_n(t) = 0$$

The transformation from one to the other is a standard technique for converting certain second-order linear differential equations into first-order nonlinear (Riccati) equations through the logarithmic derivative substitution[1].

Citations:
[1] https://en.wikipedia.org/wiki/Spherical%20Bessel%20function?useskin=monobook
[2] https://utd-ir.tdl.org/server/api/core/bitstreams/7a7995ee-a70c-4582-82e8-97cc7463b876/content
[3] https://stackoverflow.com/questions/76018812/computing-logarithmic-derivatives-of-modified-bessel-function-of-the-second-kind
[4] https://www2.compute.dtu.dk/~sohau/papers/ics2024/Computing_the_Logarithm_of_Bessel_Functions.pdf
[5] https://arxiv.org/html/2409.08729v1
[6] https://mathoverflow.net/questions/435048/second-derivative-of-the-logarithm-of-the-modified-bessel-function-of-the-first
[7] https://math.stackexchange.com/questions/2308917/does-this-inequality-about-the-second-derivative-of-logarithm-of-modified-bessel

[crowlogic]

\documentclass{article}
\usepackage{amsmath}

\title{The Logarithmic Derivative of Bessel Functions}
\author{}
\date{}

\begin{document}
\maketitle

The Bessel function $J_n(t)$ satisfies:
[
t^2J_n''(t) + tJ_n'(t) + (t^2 - n^2)J_n(t) = 0
]

Its logarithmic derivative is:
[
w(t) = \frac{J_n'(t)}{J_n(t)}
]

The logarithmic derivative satisfies a Riccati equation, as follows:

\begin{enumerate}
\item The definition directly implies:
[J_n'(t) = w(t)J_n(t)]

\item The derivative of this equation is:
[J_n''(t) = w'(t)J_n(t) + w(t)^2J_n(t)]

\item These expressions in the Bessel equation yield:
[t^2(w'(t)J_n(t) + w(t)^2J_n(t)) + t(w(t)J_n(t)) + (t^2 - n^2)J_n(t) = 0]

\item Division by $J_n(t)$ produces:
[t^2w'(t) + t^2w(t)^2 + tw(t) + (t^2 - n^2) = 0]

\item The equation rearranges to:
[w'(t) = -w(t)^2 - \frac{w(t)}{t} - (1 - \frac{n^2}{t^2})]
\end{enumerate}

This Riccati equation has the standard form:
[w'(t) = p(t) + q(t)w(t) + r(t)w(t)^2]

where:
[p(t) = -(1 - \frac{n^2}{t^2})]
[q(t) = -\frac{1}{t}]
[r(t) = -1]

\end{document}