Add support for Conversion of Ultraspherical(1) to Chebyshev()

[johannesjmeyer]

I am currently developing some code that piggybacks onto ApproxFun for the handling of spaces. Concretely, I implement a MatrixFun object that behaves like a ApproxFun.Fun, but that has matrix-valued coefficients. As I have to code up the handling of spaces myself in this case, I came across the following issue:

When I want to differentiate a (matrix-valued) Chebyshev function, I can do so by applying the matrix associated to Derivative(Chebyshev()). This, however, takes me to the space Ultraspherical(1). When I try to undo this by using Conversion(Ultraspherical(1), Chebyshev()), I get the error ArgumentError: please implement Ultraspherical(1) → Chebyshev().

Now, this is probably a bug, because
(a) ApproxFunOrthogonalPolynomials handles this case when computing derivatives itself and
(b) ApproxFunOrthogonalPolynomials.ConcreteConversion(Ultraspherical(1), Chebyshev()) actually exists.

I can use Fact (b) to construct a workaround, but I would much appreciate to keep the code clean. This should also be an easy fix on your side.

[jishnub]

I don't think ApproxFunOrthogonalPolynomials.ConcreteConversion(Ultraspherical(1), Chebyshev()) returns a meaningful answer. It's just that there aren't adequate checks to ensure that the call doesn't go through. You may check this by multiplying the conversion operators:

julia> CtoU1 = Conversion(Chebyshev(), Ultraspherical(1));

julia> U1toC = ApproxFunOrthogonalPolynomials.ConcreteConversion(Ultraspherical(1), Chebyshev());

julia> U1toC * CtoU1
ConversionWrapper : Chebyshev() → Chebyshev()
 1.0  0.0   -0.75   0.0    0.25    ⋅      ⋅      ⋅      ⋅      ⋅    ⋅
  ⋅   0.25   0.0   -0.5    0.0    0.25    ⋅      ⋅      ⋅      ⋅    ⋅
  ⋅    ⋅     0.25   0.0   -0.5    0.0    0.25    ⋅      ⋅      ⋅    ⋅
  ⋅    ⋅      ⋅     0.25   0.0   -0.5    0.0    0.25    ⋅      ⋅    ⋅
  ⋅    ⋅      ⋅      ⋅     0.25   0.0   -0.5    0.0    0.25    ⋅    ⋅
  ⋅    ⋅      ⋅      ⋅      ⋅     0.25   0.0   -0.5    0.0    0.25  ⋅
  ⋅    ⋅      ⋅      ⋅      ⋅      ⋅     0.25   0.0   -0.5    0.0   ⋱
  ⋅    ⋅      ⋅      ⋅      ⋅      ⋅      ⋅     0.25   0.0   -0.5   ⋱
  ⋅    ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅     0.25   0.0   ⋱
  ⋅    ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅     0.25  ⋱
  ⋅    ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅    ⋱

This should be the identity matrix if the operators are the inverses of each other.

Here is the actual way to obtain the conversion matrix:

julia> C = Conversion(Chebyshev(), Ultraspherical(1))
ConcreteConversion : Chebyshev() → Ultraspherical(1)
 1.0  0.0  -0.5    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   ⋅
  ⋅   0.5   0.0  -0.5    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   ⋅
  ⋅    ⋅    0.5   0.0  -0.5    ⋅     ⋅     ⋅     ⋅     ⋅   ⋅
  ⋅    ⋅     ⋅    0.5   0.0  -0.5    ⋅     ⋅     ⋅     ⋅   ⋅
  ⋅    ⋅     ⋅     ⋅    0.5   0.0  -0.5    ⋅     ⋅     ⋅   ⋅
  ⋅    ⋅     ⋅     ⋅     ⋅    0.5   0.0  -0.5    ⋅     ⋅   ⋅
  ⋅    ⋅     ⋅     ⋅     ⋅     ⋅    0.5   0.0  -0.5    ⋅   ⋅
  ⋅    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    0.5   0.0  -0.5  ⋅
  ⋅    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    0.5   0.0  ⋱
  ⋅    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    0.5  ⋱
  ⋅    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   ⋱

julia> Cinv = ApproxFun.PartialInverseOperator(C)
PartialInverseOperator : Ultraspherical(1) → Chebyshev()
 1.0  0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0  0.0  ⋯
  ⋅   2.0  0.0  2.0  0.0  2.0  0.0  2.0  0.0  2.0  ⋱
  ⋅    ⋅   2.0  0.0  2.0  0.0  2.0  0.0  2.0  0.0  ⋱
  ⋅    ⋅    ⋅   2.0  0.0  2.0  0.0  2.0  0.0  2.0  ⋱
  ⋅    ⋅    ⋅    ⋅   2.0  0.0  2.0  0.0  2.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅   2.0  0.0  2.0  0.0  2.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0  0.0  2.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0  0.0  2.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0  ⋱
  ⋮    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱   ⋱

julia> Cinv * C
TimesOperator : Chebyshev() → Chebyshev()
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
  ⋅   1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅   1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅   1.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅   1.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  ⋱
  ⋮    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱   ⋱

The inverse of a banded matrix is almost always dense, which is why the conversion is typically not defined.

Currently, there is a way to compute the derivatives in the Chebyshev basis, though:

julia> f = Fun(x->x^2, Chebyshev())
Fun(Chebyshev(), [0.5, 0.0, 0.5])

julia> f'
Fun(Chebyshev(), [0.0, 2.0])

[johannesjmeyer]

Thanks a lot for the quick reply! I will then use the partial inverse for now.

But your last example shows that the logic exists somehwere in the library. I was just wondering how to best access that logic -- as I can't just define a Fun and differentiate it, as my coefficients are matrix-valued.

[jishnub]

The last example is essentially carrying out a dense matrix-vector multiplication. This is equivalent to using the partial inverse. Maybe we can create a simpler API for such dense inverses.

[johannesjmeyer]

The dense multiplication is no issue for me at the moment, so this solution is actually quite nice. Thank you so much!