Adjoint of a polynomial not supported

[andreasvarga]

The following example shows that the transpose of a polynomial is defined, while the adjoint not,

julia> transpose(Polynomial(im))
Polynomial(im)

julia> adjoint(Polynomial(im))
ERROR: MethodError: no method matching adjoint(::Polynomial{Complex{Bool}})
Closest candidates are:
  adjoint(::Missing) at missing.jl:100
  adjoint(::Number) at number.jl:193
  adjoint(::Adjoint) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\adjtrans.jl:152
  ...
Stacktrace:
 [1] top-level scope at REPL[42]:1

I would expect the result should be Polynomial(-im).

The conjugate also fails

julia> conj(Polynomial(im))
ERROR: InexactError: Bool(-1)

[jverzani]

Thanks again. I think at one time we had the adjoint defined to be the derivative, but that is now gone and in the interim the adjoint became recursive.

I would think passing the adjoint down to the coefficients would be the right defn, right? Something like:

LinearAlgebra.adjoint(p::P)  where {P<:AbstractPolynomial} = P(adjoint(coeffs(p)), p.var)

We have conj similarly defined already:

LinearAlgebra.conj(p::P) where {P <: AbstractPolynomial} = P(conj(coeffs(p)))

But that fails as the type of P is Complex{Bool} but needs to be Complex{Int64}. Will make the necessary adjustment (along with passing the variable symbol through).

Anyways, there is a choice between adjoint.(coeffs(p)) versus adjoint(coeffs(p)). I think the latter is the correct one.

[andreasvarga]

I discovered the above issues by manipulating polynomial matrices. By the way, I registered my recent package MatrixPencils.jl , where the Polynomial matrix type is fully supported as an alternative to 3-d array representation of polynomial matrices. I wonder if at certain moment, it could be an issue for you to support some basic operations for polynomial matrices as well. For such a matrix, the variable could be defined globally for all entries.

For me the next step will be to handle in my software rational matrices. These are matrices whose entries are rational functions (i.e., quotients of two polynomials). For the beginning I will use two 3-d matrices to store the numerator and denominator polynomials. Do you plan any extension of the polynomial type to rational function type (e.g., using a touple of two polynomials, representing the numerator and denominator)? This would be a great alternative way to handle rational matrices.

[jverzani]

I'd love to make Polynomials as generic as possible, so supporting matrices is something to look at. If I have it right, the first issue would be relaxing the AbstractPolynomial{T <: Number} restriction, which is actually pretty baked in, so this would be an undertaking to see if it could be done. (I actually removed the T<:Number restriction in the base type, but it is still there for subtypes.)

As for rational function types. I think that might be best in a separate package. (Then also the p ade functionality could get some needed attention.) Starting with a tuple or struct for the numerator and denominator seems like the right starting point.

[jverzani]

Closed by #226

[hurak]

May I still add my viewpoint on the issue of conj() for polynomials?

At the end of the day, what is actually (complex-)conjugated is the VALUE of the polynomial at some "point". If the polynomial is to be evaluated on the real line, then the operation that turns the original polynomial into one that evaluates to the conjugate value is indeed the conjugation of the coefficients.

For

a(x) = a₀ + a₁ x + a₂ x²

the conjugate of the value of the polynomial evaluated at 3 is

conj(a(3)) = conj(a₀) + conj(a₁) 3 + conj(a₂) 3²

Therefore the conjugate of the polynomial is given

̃a(x) = conj(a₀) + conj(a₁) x + conj(a₂) x² 

However, if the polynomial is to be evaluated on the imaginary axis of the complex plane, something else must be done with the original polynomial so that it evaluates to the conjugate of the value of the original polynomial.

If the polynomial a is evaluated at 3im

conj(a(3im)) = conj(a₀) + conj(a₁) conj(3im) + conj(a₂) conj(3²im²) = conj(a₀) - conj(a₁) conj(3im) + conj(a₂) (3²im²)

Therefore the conjugate of the polynomial is given

̃a(s) = conj(a₀) - conj(a₁) s + conj(a₂) s² 

Let me simplify this for polynomials with real coefficients

̃a(s) = a₀ - a₁ s + a₂ s² 

The pattern is obvious: the coefficients with odd powers of the variable change sign.

Compactly (for a real-coefficient polynomial),

̃a(s) = a(-s) 

In the above, I intentionally switched the name of the variable from x to s because such is the convention in the fields of signal processing, filter design and control theory (possibly p too), in which (fractions of) polynomials are used to describe continuous-time linear filters and controllers. Behaviour of these on the imaginary axis is often analyzed.

Finally, let me add into the confusion. In the just mentioned domains of signals and (control) systems, there might also be a need to evaluate the polynomials on the unit circle in the complex plane. In fact, my guess is that this is the practically most useful situation of all the three just enumerated (real axis, imaginary axis, unit circle). It is very much relevant in design of discrete-time (digital) filters and controllers. For a polynomial

a(z) = a₀ + a₁ z + a₂ z²

in which the variable is now z (related to z-transform), the conjugate of the value of the polynomial evaluated at some particular z=exp(im ω) is

conj(a(exp(im ω))) = conj(a₀) + conj(a₁) exp(-im ω) + conj(a₂) exp(-2 im ω)

conj(a(z)) = conj(a₀) + conj(a₁) 1/z + conj(a₂) 1/z²

Compactly (again, for real coefficients)

̃a(z) = a(1/z) 

This can be rewritten as

̃a(z) =a₀ + a₁ 1/z + a₂ 1/z² 

In some domains this is actually called paraconjugate, see for example https://ccrma.stanford.edu/~jos/filters/MIMO_Paraconjugate.html

A funny thing is that even though the resulting is a rational function and not a polynomial, it can be viewed as a polynomial too by invoking the concept of a Laurent polynomial. And we now have it in Polynomials package :-) The variable of the (para)conjugate polynomial is then 1/z. I am still not quite fluent with all the Julia tricks and rules but perhaps the unit-circle version of the conjugation (paraconjugation) should only be defined for Laurent polynomials to guarantee type stability.

To summarize, different occassions might call for different versions of conjugate of a polynomial. Now what? One possibility is to give them different names. The third version (the one on the unit circle) seems to have one (adopted by at least some community). Another possibility to choose the version of the conjugate operation base on the name of the variables. Say, by default the real variable version is considered, and if the variable is s or p, the imaginary variable is considered, and if the variable is z or λ (or perhaps q and d, which are also common), the unit-circle version is used.

Just a suggestion.

[andreasvarga]

When I raised the issue withconj, I already had a bad feeling, because I was aware of the involved possible complications. So, many thanks for addressing these issues. In symbolic languages, the standard usage of conj for a polynomial is simply to take the conjugates of coefficients and replace the variable, say s, by conj(s). For example,

P = s^3 - 3*s^2 + 3*s - 1

conj(P) = conj(s)^3 - 3*conj(s)^2 + 3*conj(s)  - 1

I see three possibilities:

1. Julia allows to add an optional argument for conj (or even a keyword argument), by specifiying the involved operation on the variable (e.g., s -> -s or z -> 1/z). I think this is the most flexible and general solution and could also cover more general Moebius transformations (provided the rational function type will be sometime supported as well).

2. Using special variable names ass (the Laplace transform variable) orz (the Z-transform variable) could be an alternative for setting default options. This option could occasionally lead to misleading results, since we can not assume that all users are aware of these assumed meanings of the variable names.

3. Using different names asconj(to conjugate only the coefficients, but not the variable), cconj (to conjugate the coefficients and perform s -> -s), and dconj (to conjugate the coefficients and perform z -> 1/z). This option is probably the easiest to implement.

[hurak]

And one more thing. Perhaps it is wise to insert it into this issue instead of opening a new one, because it can be viewed as an extension of what has been discussed here - matrix polynomials.

By a matrix polynomial we understand a univariate polynomial whose coefficients are matrices, such as

        [ 1  1 ]   [ 2  0 ]     [ 3  0 ]
P(s) =  |      | + |      | s + |      | s²
        [ 0  0 ]   [ 2  1 ]     [ 0  0 ].

A very interesting (in fact, rather fundamental) observation is that such object can also be viewed as a polynomial matrix, that is, a matrix whose entries are univariate polynomials

        [ 3s²+2s+1     1 ]
P(s) =  |                |
        [    2s        s ].

Let me inform you that there is a Julia package for polynomial matrices called PolynomialMatrices.jl. It was developed by Niklas Everitt @neveritt quite some time ago (around 2016-2017). I stepped in a year ago or so, trying to save the then essentially broken package by upgrading it to Julia v1. As of now, apparently, the package has not yet incorporated the breaking change(s) in v0.8 Polynomials package (but this could perhaps be easy and I will do it soon.)

The package aims at providing some basic operations for matrices of univariate polynomials such as determinant and various decompositions and triangularizations. Mostly motivated by the signals and (control) systems fields. The motivation was/is to bring into Julia the functionality provided by the commercial Polynomial Toolbox for Matlab.

The thing is, however, that even if the resulting object is viewed as a matrix (of polynomials), it is conveniently entered as a (matrix) polynomial. In fact, this is also how the polynomial matrix is represented internally in PolynomialMatrices. See the README on the front page, where we explain it in detail. Technically speaking, the matrix is stored internally as a dictionary. This would correspond to SparsePolynomial type in your Polynomials package, I guess. But entering the polynomial in this way (see the README) is quite akward, I am afraid (in fact, it was @andreasvarga question/request some time ago that made me to realize it). An alternative way for entering polynomials provided by the package is through a 3D array, but this is not much more user friendly. What would be absolutely perfect is to be able to do something like this

using Polynomials

A0 = rand(2,2)
A1 = rand(2,2)
A2 = rand(2,2)

s = Polynomial([0,1],:s)

A = A0 + A1*s + A2*s^2

Whether it is wise to implement this in the core Polynomials package (it could make sense because the s variable would come from) or whether it deserves its own package (perhaps the PolynomialMatrices package after some upgrade), I do not know. I have not mastered Julia enough to come up with a true Julian design decision by myself. That is why I am presenting the issue here.

[jverzani]

Thanks for the attention here. I'd have to look up the details, but there was some discussion about how to define conj before we settled on what we have done. I personally was looking for something with a reasonably well-known expectation. With the LaurentPolynomial type, we could certainly provide an alternate definition for conj. I think this makes sense at first glance (it is available, easy to document for users, has a mathematical reasoning, ...). I would be really wary of some keyword argument or special convention regarding symbols. If people wanted such for the standard basis polynomials, a new method name would be cleanest, I think.

As for matrix based coefficients. That will require some work. At one point recently I was coding up Bernstein polynomials and thought it would be cute to have vector-valued coefficients for a spline example. I could get that to work only by working around some of the generic assumption for the underlying type T in the common.jl file. So possible, but does require some effort.

Were this done, the natural constructor would be to use a vector of matrices to specify the coefficients ([ rand(2,2), rand(2,2), rand(2,2) ]). I'm not grasping the need for a 3d array here. The SparsePolynomial class could be used to specify this with a dictionary of matrices (e.g. Dict(0=>rand(2,2), 1=>rand(2,2), 2=>rand(2,2)) and things should flow from there for polynomials over the set of matrices.

[hurak]

Were this done, the natural constructor would be to use a vector of matrices to specify the coefficients ([ rand(2,2), rand(2,2), rand(2,2) ]). I'm not grasping the need for a 3d array here.

Indeed. I think that the current implementation in PolynomialMatrices, which works like this

julia> A = zeros(Int,2,2,3);
julia> A[:,:,1] = [1 1;0 0];
julia> A[:,:,2] = [2 0;2 1];
julia> A[:,:,3] = [3 0;0 0];

julia> P = PolyMatrix(A,:s)
2×2 PolyArray{Int64,2}:
  Poly(1 + 2⋅s + 3⋅s^2)  Poly(1)
  Poly(2⋅s)              Poly(s)

may just reflect the Matlab style of thinking about arrays of matrices. But a more Julian way would be, of course, to use a vector of matrices.

The SparsePolynomial class could be used to specify this with a dictionary of matrices (e.g. Dict(0=>rand(2,2), 1=>rand(2,2), 2=>rand(2,2)) and things should flow from there for polynomials over the set of matrices.

This is how PolynomialMatrices have it right now:

julia> d = Dict(0=>[1 1;0 0], 1=>[2 0;2 1], 2=>[3 0;0 0]);
julia> P = PolyMatrix(d,:s)
2×2 PolyArray{Int64,2}:
  Poly(1 + 2⋅s + 3⋅s^2)  Poly(1)
  Poly(2⋅s)              Poly(s)

[jverzani]

I opened an issue #231 and PR #232 which attempts to loosen the T <: Number assumption in the Polynomial type. It is only partially successful.

[hurak]

Let me inform you here that I started a small project aiming at implementing some solvers for some equations with polynomials: PolynomialEquations.jl. In particular, you may want to have a look at solvers for symmetric conjugate equations, which are a perfect example of a situation where conjugation of a polynomial with respect to the imaginary axis is needed. The unit-disc version is in to-do list. I have implemented both cconj() and dconj() functions for this purpose.

[jverzani]

Thanks, let me have a look.

I kinda put this down for a bit after working on the PR to allow polynomials with matrix coefficients. I'd be curious if you thought that was helpful at all.

[neveritt]

Thanks for the nice discussions here. I would just like add to the possible implementations of conj mentioned by @andreasvarga. A fourth possibility is explored in RationalFunctions.jl where a type parameter is used to keep track if the variable is conjugated or not.

As @jverzani mentioned the earlier discussion on conj for polynomials can be found here

[jverzani]

In #242 I added paraconj and cconj for use with the LaurentPolynomial type. I did not add the fourth possibility mentioned here, though if pushed it could be.