Nth Derivative, but N is a fraction

Half Derivatives: An Operator Theory Perspective

Intro

When we’re taught calculus, we’re taught second derivatives are a derivative of a derivative, and the third derivative is the derivative of that and so on.

A natural question is: Is there a “half derivative;” an operator which if you apply it twice it amounts to applying a derivative just once? Concretely, is there an operator H for which H(Hf) = Df? The answer is yes, at least for smooth functions f.

Furthermore, you can define a 1/7 derivative (which if you apply seven times, you get a first derivative) or a 2/3 derivative (which if you apply three times, you get a second derivative). In fact for any rational number, m/n, you can define the m/n derivative with the property that:

I’ll call this property the “implicit” property. Operators that satisfy this implicit property are called “fractional derivatives.” Throughout the paper, we will at times focus specifically on the half-derivative, but at other times will point out that the calculation can be done with more general fractions.

Actually this can be defined for irrational numbers as well, via limits of rational numbers. Fractional derivatives are weird, but they do have a nice property that they interpolate the whole number derivatives that we already know and love. In this picture the green line goes continuously from the 0th derivative of 1/2*x² to its 2nd derivative.

Image source: Wikimedia Commons.

This blog will walk through the calculation of the fractional derivative of some common functions:

• The exponential

• Sine and cosine

• Polynomials

Though there are lots of ways to talk about this topic, I’ve chosen to structure this to touch on a few math topics, making potentially surprising connections.¹

Recommended Background

Assumes some familiarity with calculus and linear algebra.

Calculating the Half Derivative: Special Functions

In general half derivatives can be difficult to calculate, but let’s look at one case where it’s not difficult: The exponential function is an “eigenfunction” of the derivative, meaning that its derivative is the same function up to a constant multiplier.

\(D^n\left(e^{\lambda x}\right)=\lambda^ne^{\lambda x}\)

You may then guess that we can replace the whole number n with any real number α.

\(D^\alpha\left(e^{\lambda x}\right)=\lambda^\alpha e^{\lambda x}\)

If you check, you’ll see that an operator defined in this way satisfies the implicit property for the fractional derivative. For example:

\(D^{2/3}D^{2/3}D^{2/3}\left(e^{\lambda x}\right)= D^{2/3}D^{2/3}\left(\lambda^{2/3}e^{\lambda x}\right)= D^{2/3}\left(\lambda^{4/3}e^{\lambda x}\right)= \lambda^2e^{\lambda x} = D^2\left(e^{\lambda x}\right)\)

But what about a function like sine? An expert trigonometer may remember a helpful trig identity here:

\(\sin(x) = -\frac{i}{2}e^{ix}+\frac{1}{2}e^{-ix}\)

Knowing then what we know about the fractional derivative of the exponential, we can conclude that:

\(D^\alpha\sin(x)=-\frac{i}{2}i^\alpha e^{ix}+\frac{i}{2}\left(-i\right)^\alpha e^{-ix}\)

But check this out: Using the identity that i=e^(iπ/2) and -i=e^(-iπ/2), we can write this as:

\(\begin{array}{rcl} D^\alpha\sin(x) &=& -\frac{i}{2}e^{i\alpha\pi/2}e^{ix} +\frac{i}{2}e^{-i\alpha\pi/2}e^{-ix} \\ &=& -\frac{i}{2}e^{i(x+\alpha\pi/2)} +\frac{i}{2}e^{-i(x+\alpha\pi/2)} \\ &=& \sin(x+\alpha\pi/2) \end{array}\)

This is a really cool representation, because we see that at α=1, this is cos(x), and as α increases from 0 to 1, the graph slowly changes phase.

We got this fractional derivative by writing sine as a sum of exponentials. Cosine can be handled similarly. In fact this representation shows that sine and cosine are just pit stops on a pair of continuous rotations through the complex unit circle. For other functions, we aren’t as lucky. Most functions can’t be written as a sum of exponentials, but this example motivates our approach.

The last function that we promised to calculate a half-derivative of is the power function x^m. But to do this will take more work. In particular, we will need to take the Fourier transform. If you want to take this approach as a lucky guess, you can jump to the last section. Otherwise we will take a long detour in an attempt to explain why the Fourier transform may be a natural thing to try.

Matrix Operators

Before we go forward, let’s first talk about matrices. The reason is that we’ve been talking about derivatives as linear operators, which makes them generalizations of matrices. (We’ve seen this analogy already when we mentioned eigenfunctions.) Basic operator theory is motivated from things known about matrices. For example, the half derivative is the square root of the derivative in the matrix sense: If the square root of a matrix M is A then A(Av) = Mv for any vector v. Similarly if we call the square root of the derivative “H”, then H(Hf) = Df; so H is the half-derivative.

---

Before we go on, I want to introduce bra-ket notation for matrices. We write a vertical bar followed by angle bracket for a vector, and an angle bracket followed by a vertical bar for it’s conjugate transpose.

\(|\varphi\rangle=\left(\begin{array}{c}a\\b\end{array}\right),\quad\quad\langle\varphi|=\left(\begin{array}{cc}a&b\end{array}\right)\)

This makes writing inner and outer products easy, for example.

\(\langle\varphi|\varphi\rangle=a^2+b^2,\quad\quad |\varphi\rangle\langle\varphi|=\left(\begin{array}{cc} a^2 & ab \\ ab & b^2 \end{array}\right)\)

---

So how exactly do we take the square root of a matrix? One way is with the Taylor series.

Unfortunately the Taylor series for the square root doesn’t exist, because √x isn’t differentiable at zero. We use a common trick, where we instead think about the Taylor series of √(1+x). We will worry about the “1+” later, but for now, we recall/lookup the Taylor series.

\(\sqrt{1+x} = 1+\frac{1}{2}x-\frac{1}{8}x^2+\frac{1}{16}x^3-\frac{5}{128}x^4+\cdots\)

This Taylor series is defined for real numbers, x. But it ought to also work for matrices; that is, for a matrix, M, it ought to hold that:

\(\sqrt{I+M} = I+\frac{1}{2}M-\frac{1}{8}M^2+\frac{1}{16}M^3-\frac{5}{128}M^4+\cdots\)

Now √(I+M) should be a matrix which, if multiplied by itself, returns I+M. It’s not obvious that the Taylor series, which worked for number variables, will work the same way for matrices, but it’s a good guess.

To verify that it’s a correct guess, we can take the Taylor series and multiple it by itself. The resulting Taylor series should equal I+M, which would tell us that we indeed have the square root. To multiple two infinite series, we have to do a lot of distributing. By writing in a grid and multiplying all terms by all terms, we have that all terms containing a certain power of M lie in the same diagonal. For example, in the grid below all the powers of M³ fall in the purple diagonal. (We’ve only filled out the grid up to powers of of 4.)

We see that the first two diagonals add to I (the yellow diagonal) and M (the green diagonal). The remaining diagonals cancel out and are all zero! This is a fun calculation to see, but actually not surprising because this is what would happen if we did this multiplication on real number Taylor series.

A similar Taylor series could be done for matrices raised to other powers.

Singular Value Decomposition

If you knew already that Taylor series apply to matrices, then this is a natural first approach to taking a square root. But Taylor series are difficult to caclulate; in practice you could only ever get an approximation, because you can’t calculate all the terms. There is however a more clever way to take the square root, which we will explore in this section.

To get there, we take the singular value decomposition (SVD). Though the SVD has a more general form than we cover, we will talk about a simplified version that only applies to Hermitian matrices; these are matrices that are equal to their own conjugate transpose.

The SVD says: If a Hermitian matrix, M, has unit eigenvectors φ_i with corresponding eigenvalues λ_i, then it can be written as:

\(M = \sum_i\lambda_i|\varphi_i\rangle\langle\varphi_i|\)

Another similar formula (not part of the SVD) that will be useful later is:

\(I=\sum_i|\varphi_i\rangle\langle\varphi_i|\)

This is getting very abstract at this point. So let’s focus on a specific matrix.

\(\left(\begin{array}{cc} 5&4\\4&5 \end{array}\right)\)

This matrix has eigenvectors φ_1 = (1/√2, -1/√2) and φ_2 = (1/√2, 1/√2), with eigenvalues λ_1 = 1 and λ_2 = 9. You can check this. And you can verify that therefore,

\(\left(\begin{array}{cc} 5&4\\4&5 \end{array}\right)=1|\varphi_1\rangle\langle\varphi_1|+9|\varphi_2\rangle\langle\varphi_2|\)

where the outer products are

\(|\varphi_1\rangle\langle\varphi_1|=\left(\begin{array}{cc} 1/2&-1/2\\-1/2&1/2\end{array}\right)\quad\text{and}\quad |\varphi_2\rangle\langle\varphi_2|=\left(\begin{array}{cc} 1/2&1/2\\1/2&1/2\end{array}\right)\)

But watch what happens when we square the equation.

\(\begin{array}{rcl} \left(\begin{array}{cc}5&4\\4&5\end{array}\right)^2&=&\left(1|\varphi_1\rangle\langle\varphi_1|+9|\varphi_2\rangle\langle\varphi_2|\right)^2\\ &=&1^2|\varphi_1\rangle\langle\varphi_1|\varphi_1\rangle\langle\varphi_1| +\\ &=&1\cdot9|\varphi_1\rangle\langle\varphi_1|\varphi_2\rangle\langle\varphi_2| +\\ &=&9\cdot1|\varphi_2\rangle\langle\varphi_2|\varphi_1\rangle\langle\varphi_1| +\\ &=&9^2|\varphi_2\rangle\langle\varphi_2|\varphi_2\rangle\langle\varphi_2| \end{array}\)

If you look the middle of each term ends up being an inner product! But φ_1 and φ_2 are orthonormal: The inner product of φ_1 and φ_2 is 0, but the inner product of φ_i and φ_i is 1. So the equation simplifies dramatically to:

\(\left(\begin{array}{cc}5&4\\4&5\end{array}\right)^2=1^2|\varphi_1\rangle\langle\varphi_1|+9^2|\varphi_2\rangle\langle\varphi_2|\)

The square just got applied to the eigenvalues! You can again verify this is true.

\(\left(\begin{array}{cc}5&4\\4&5\end{array}\right)^2=1\left(\begin{array}{cc}1/2&-1/2\\-1/2&1/2\end{array}\right)+81\left(\begin{array}{cc}1/2&1/2\\1/2&1/2\end{array}\right)=\left(\begin{array}{cc}41&40\\40&41\end{array}\right)\)

But you can see if you multiply by the matrix again, the same thing will happen again, so you generally have:

\(\left(\begin{array}{cc}5&4\\4&5\end{array}\right)^n=1^n|\varphi_1\rangle\langle\varphi_1|+9^n|\varphi_2\rangle\langle\varphi_2|\)

For this remarkable property, all we needed was for the vectors making up the outer products to be orthonormal. And eigenvectors are always orthonormal. So we have a rule for any matrix with its eigenvalues and eigenvectors.

\(M^n=\sum_i\lambda_i^n|\varphi_i\rangle\langle\varphi_i|\)

You may see where this is going: If I can do this for any whole number n, then I ought to be able to do it for any real number α. And that’s absolutely right. But let’s check why that’s true. Watch these equations below: We just rewrite √(I+M) with the Taylor series, then we can replace each M^n with the same equation on λ^n. After we factor out the outer products, this becomes the same Taylor series on λ, so we can replace back with √(I+λ).

\(\begin{array}{rcl} \sqrt{I+M}&=&I+\frac{1}{2}M-\frac{1}{8}M^2+\frac{1}{16}M^3+\cdots\\ &=&\sum_i|\varphi_i\rangle\langle\varphi_i|+\frac{1}{2}\sum_i|\lambda_i\varphi_i\rangle\langle\varphi_i|-\frac{1}{8}\sum_i\lambda^2|\varphi_i\rangle\langle\varphi_i|+\cdots\\ &=&\sum_i\left(1+\frac{1}{2}\lambda_i-\frac{1}{8}\lambda_i^2+\frac{1}{16}\lambda_i^3+\cdots\right)|\varphi_i\rangle\langle\varphi_i|\\ &=&\sum_i\sqrt{1+\lambda_i}|\varphi_i\rangle\langle\varphi_i| \end{array}\)

I’ll now tell you that you can multiply I in that equation by some small number, ε, and get the same result.

\(\sqrt{\varepsilon I+M}=\sum_i\sqrt{\varepsilon+\lambda_i}|\varphi_i\rangle\langle\varphi_i|\)

If your careful with the analysis, you can let that small number approach zero, and you get the desired result that:

\(\sqrt{M}=\sum_i\sqrt{\lambda_i}|\varphi_i\rangle\langle\varphi_i|\)

Checking on our matrix from above, we get:

\(\left(\begin{array}{cc}5&4\\4&5\end{array}\right)^{1/2}=1\left(\begin{array}{cc}1/2&-1/2\\-1/2&1/2\end{array}\right)+3\left(\begin{array}{cc}1/2&1/2\\1/2&1/2\end{array}\right)=\left(\begin{array}{cc}2&1\\1&2\end{array}\right)\)

And calculating that was much easier than calculating the Taylor series! Again, you can do the same trick with any power of a matrix, and again get the power to apply to the eigenvalue coefficients.

\(M^\alpha=\sum_i\lambda_i^\alpha|\varphi_i\rangle\langle\varphi_i|\)

In fact, any function with a power series, can be written this way.

Fourier Transform

Returning now to the derivative. Derivatives are linear operators on the vector space of smooth functions. But unlike 2-D vectors like we used in our example, function spaces are infinite dimensional, meaning that they don’t have any finite basis.

The SVD on matrices is a summation over a finite number of eigenvectors. In this section we make an analogy with an infinite number of eigenfunctions.

For matrices, we worked with eigenvectors as a basis. For operators, let’s choose a basis composed of eigenfunctions. Let’s try taking this family of eigenfunctions: {exp(2πiξx) : ξ a real number}.

The first question you should ask is: Can you write functions as a combination of this family of functions. For matrices, we wrote functions as summations of basis vectors. But now, we need to write an integral, because there’s a continuous spectrum of basis functions. We write the coefficients as a function, f-hat, to emphasize that there’s a coefficient for every real number. In the right-most expression, we write exp(2πiξx) as a vector to remind that it lives in a vector space.

\(f(x)=\int_{-\infty}^\infty\hat f(\xi)e^{2\pi ix\xi}d\xi=\int\hat f(\xi)\left|e^{2\pi i x\xi}\right\rangle\)

This may look familiar; f-hat is the Fourier transform! It turns out f can be written this way if it is integrable; that is, if it has a finite integral over all real numbers.

The second question you should ask is: Are these eigenfunctions orthonormal? Recall that the standard inner product on functions is defined as:

\(\langle f|g\rangle=\int_{-\infty}^\infty f^*(x)g(x)dx\)

where the * means the complex conjugate.

The answer here is: Yes, but… Rather than get into the gory details of that “but,” let’s instead state the Fourier inversion theorem, which says that:

\(\hat f(\xi)=\int_{-\infty}^\infty f(x)e^{-2\pi i\xi x}dx\)

Another way to say this is:

\(\hat f(\xi)=\left\langle e^{2\pi i\xi x}|f\right\rangle\)

The sign change in the exponent is because the left side version of the bra-ket notation does the complex conjugate. What this is saying though is that applying

\(\left\langle e^{2\pi i\xi x}\right|\)

to the left side of our equation above “cancels” all the terms except for the ξ-th term. This is exactly what we desire from orthonormality. So this is enough

Above we called f-hat the Fourier transform. Formally the Fourier transform, marked with a calligraphic F, is the operator that takes f to f-hat.

\(\mathcal F(f):=\hat f\)

Let’s see what happens when we apply the Fourier transform on the derivative.

\(\mathcal F(Df)=\int_{-\infty}^\infty Df(x)e^{-2\pi i\xi x}dx\)

We’re going to do integration by parts with the function f and exp(-2πiξx).

\(\mathcal F(Df)=\int_{-\infty}^\infty Df(x)e^{-2\pi i\xi x}dx =\left. f(x)e^{-2\pi i\xi x}\right|_{x=-\infty}^\infty +2\pi i\xi\int_{-\infty}^\infty f(x)e^{-2\pi i\xi x}dx\)

The first term must approach zero for both + and -infinity. This isn’t obvious, but it turns out that if this isn’t true, then f isn’t going to zero fast enough to be integrable. So this becomes simply:

\(\mathcal F(Df)=2\pi i\xi \mathcal F(f)\)

Before we extended this, let’s notice what is saying. Because Df is just a function, we can write it as:

\(Df(x)=\int_{-\infty}^\infty\widehat{Df}(\xi)e^{2\pi ix\xi}d\xi\)

But what we just showed is that the coefficients are just simply 2πiξ times the coefficients we’d have if we just put in f.

\(D|f\rangle = \int 2\pi i\xi\hat f(\xi)\left|e^{2\pi ix\xi}\right\rangle=\int 2\pi i\xi\left|e^{2\pi ix\xi}\right\rangle\hat f(\xi)\)

But remember from above that f-hat is just a particular inner product on f.

\(D|f\rangle=\int 2\pi i\xi\left|e^{2\pi ix\xi}\right\rangle\left\langle e^{2pi ix\xi}|f\right\rangle\)

To write D as an operator, we can just leave the f off.

\(D=\int 2\pi i\xi\left|e^{2\pi ix\xi}\right\rangle\left\langle e^{2pi ix\xi}\right|\)

Written like this, we now see how this generalizes the SVD. Not surprisingly the coefficient is the corresponding eigenvalue to the eigenvector.

Continuing, simple repeated integration by parts tells us that:

\(\mathcal F(D^nf)(\xi) = \left(2\pi i\xi\right)^n\mathcal F(f)(\xi)\)

We can do the usual trick to extend to other powers. For square root, for instance, we have:

\(\begin{array}{rcl} \mathcal F\left(\sqrt{I+D}\right)(f)&=&\mathcal F(f)+\frac{1}{2}\mathcal F(Df)-\frac{1}{8}\mathcal F(D^2f)+\frac{1}{16}\mathcal F(D^3f)+\cdots\\ &=&\mathcal F(f) + \frac{1}{2}\left(2\pi i\xi\right)\mathcal F(f) -\frac{1}{8}\left(2\pi i\xi\right)^2\mathcal F(f) + \frac{1}{16}\left(2\pi i\xi\right)^3\mathcal F(f) +\cdots\\ &=&\sqrt{1+2\pi i\xi}\mathcal F(f) \end{array}\)

Using the appropriate Taylor series for the power, then, you can conclude:

\(\mathcal F\left(D^\alpha f\right)(\xi)=\left(2\pi i\xi\right)^\alpha\mathcal F(f)(\xi)\)

Fractional Derivatives: Polynomials

We’re now finally ready to calculate the fractional derivative of the power function, x^m. We want to do this calculation for a general m, so we will only consider positive x, because x^m is not defined when x is negative and m is irrational.

Let’s start with the Fourier transform of f(x)=x^m, for x > 0.

\(\mathcal F(f)(\xi) = \int_0^\infty x^me^{-2\pi i\xi x}dx\)

The lower bound is zero, because we only consider positive values of x, and define f to be zero elsewhere. As an integration rule, whenever we see x appear as both a base and an exponent in a integrand, we should try a gamma function. And actually only a simple substitution is needed to get:

\(\mathcal F(x^m)(\xi)=\frac{\Gamma(m+1)}{\left(2\pi i\xi\right)^{m+1}}\)

The Fourier transform of a power function on positive x is a function with ξ only in the denominator with an exponent of at least 1. Because we can multiply the input and the output by the same constant (remember it’s just an integral), any function that transforms to have a ξ only in the denominator with an exponent of at least 1 must be a power function. We’ll take for granted that two smooth functions with the same transform are the same.

Using the rule we made in the last section, we have for any α < m:

\(\mathcal F\left(D^\alpha(x^m)\right)(\xi)=\left(2\pi i\xi\right)^\alpha\frac{\Gamma(m+1)}{\left(2\pi i\xi\right)^{m+1}} = \frac{\Gamma(m+1)}{\left(2\pi i\xi\right)^{m-\alpha+1}}\)

Because α < m, this again has a ξ only in the denominator with an exponent of at least 1. So it’s the Fourier transfer of a different power function. Specifically,

\(\mathcal F\left(D^\alpha(x^m)\right)(\xi)= \frac{\Gamma(m+1)}{\left(2\pi i\xi\right)^{m-\alpha+1}} =\frac{\Gamma(m+1)}{\Gamma(m-\alpha+1)}\mathcal F\left(x^{m-\alpha}\right)(\xi)\)

This tells us that:

\(D^\alpha(x^m)=\frac{\Gamma(m+1)}{\Gamma(m-\alpha+1)}x^{m-\alpha}\)

Look at how cool this formula is, though! If you’ve encountered the gamma function before, then you know that it interpolates the factorial. When m and α are both integers then this boils down to:

\(D^\alpha(x^m)=\frac{m!}{(m-\alpha)!}x^{m-\alpha}\)

This is formula you already know for repeated derivatives of power functions. If you had guessed that the way to generalize to fractional derivatives is to replace the factorials with the gamma function, your guess would have been exactly correct!

Other Things to Think About

To review: We found with our earlier examples (exponentiation and sine) that a fractional derivative is easy to compute when you’re able to write a function as a linear combination of eigenfunctions. A detour through matrices showed why eigenvectors are so useful for this. Though we lucked out with those two functions, we wanted a more general way to write a function as a combination of eigenfunctions. To do this we looked at the Fourier transform. And we finally used this to calculate the fractional derivative of the power function.

These techniques help explore fractional derivatives, which are interesting in their own right, but the techniques work for smooth functions (more general than powers) and for operators (more general than derivatives).

Some other things I’ve been thinking about:

• Do fractional derivatives provide a homotopy of I and D in some function space?

• We used eigenfunctions of D. What are the eigenfunctions of partial derivatives? What about fixed points?

• What is the fractional derivative on logarithms? Other functions?

To read more about fractional derivatives (more concrete properties and an application), I recommend this blog. And you can also check out the wiki.

Thanks to Melanie Beck and Sisi Song for feedback.

1. This is a copy of a blog original written here.