Some thoughts on conjugation

TL;DR

The package has some problems related to the conjugation of Poly objects. In my opinion, conj should not only return a new Poly with conjugated coefficients but also keep track of its input data. The solution proposal in #59 seems to be true only for real inputs.

What I would expect from conj is that conj(a*x) = conj(a)*conj(x). In fact, after obtaining conj of a Poly, the variable x has been also changed to y = conj(x). What I mean by that is conj(p(x)) + p(x) should not be allowed for a Poly object p(x) with the current implementation.

Read below for more...

Example

julia> p3 = Poly([1+1im, 1])
Poly((1 + 1im) + x)
julia> p4 = conj(p3)
Poly((1 - 1im) + x)

julia> p3(1)
2 + 1im
julia> p4(1)
2 - 1im
julia> conj(p3(1)) == p4(1)
true

julia> p3(1 + 1im)
2 + 2im
julia> p4(1 + 1im)
2 + 0im
julia> conj(p3(1 + 1im)) == p4(1 + 1im)
false

To me, the variable name together with its conjugation property should fix the variable information of a Poly object. Actually, I have implemented RationalFunctions.jl in this way, and that seems to be working.

The way I implemented RationalFunction objects was to encode both the variable name (in this case, p.var) and its conjugation property (for example, a value type such as Conj{true}) into the type itself. Then, I could rely on the dispatching mechanism of Julia for all mathematical operations without even caring about any if checks to guarantee same-variable-polynomial operations. I am not sure if this is the best practice, but it seems to be working.

I believe,

[x] It provides safe operations without the need for explicitly checking in run-time if the corresponding rational functions (in this case, polynomials) are of the same-variable,
[x] It allows for keeping track of conjugated variables,
[x] It allows for efficient function calls (i.e., (r::RationalFunction{T,Conj{true}}){T}(x::Real) = _funceval(r, x)) via the dispatching mechanism,
[ ] It pollutes (a bit ?!) more the source code, and maybe makes it harder to read, and,
[ ] It might get slower if one needs several different variables and conjugation,
- However, I doubt that a user will have a lot of variables. In such a case, one can also provide some precompilaton procedure for the usual suspects (such as :x, :t, :s, :q, etc.).

Anyways, in the current implementation of Polynomials, p.var is hardly ever used except for checking for consistent operations, and displaying purposes. If above functionality (i.e., tracking of conjugated variables) is to be implemented without breaking the coefficients' structure (in this case, Vector or T<:AbstractVector if sparse coefficients are going to be supported), the other viable option is to have a Bool in the type structure. Then, this will bring another if switch into the picture.

I just wanted to share my opinions about the issue, as I am interested in having RationalFunctions with real coefficients and complex inputs in my applications. Right now I seem to have solved it properly in the above-mentioned package. However, if anyone is interested in some brainstorming and/or unifying these representations so that both the packages could play well with each other, I am open to discussions.

Apart from this, I have also a couple of concerns regarding the convenient ctranspose for polyder, and the no-op transpose deprecation warning in Julia 0.5.

Example

julia> using Polynomials

julia> p1 = Poly([1,1])
Poly(1 + x)

julia> p2 = Poly([2,2])
Poly(2 + 2⋅x)

julia> A = [p1 p2; p2 p1]
2×2 Array{Polynomials.Poly{Int64},2}:
Poly(1 + x)    Poly(2 + 2⋅x)
Poly(2 + 2⋅x)  Poly(1 + x)  

julia> b = [p1, p2]
2-element Array{Polynomials.Poly{Int64},1}:
Poly(1 + x)  
Poly(2 + 2⋅x)

julia> A*b # perfectly fine
2-element Array{Polynomials.Poly{Int64},1}:
Poly(5 + 10⋅x + 5⋅x^2)
Poly(4 + 8⋅x + 4⋅x^2) 

julia> b.'*A*b # Warning on no-op `transpose`
1-element Array{Polynomials.Poly{Int64},1}:
Poly(13 + 39⋅x + 39⋅x^2 + 13⋅x^3)

julia> b'*A*b # Different result from above (of course!)
1-element Array{Polynomials.Poly{Int64},1}:
Poly(13 + 26⋅x + 13⋅x^2)

julia> b'*b # MethodError for `dot(p1::Poly, p2::Poly)`
ERROR: MethodError: no method matching dot(::Polynomials.Poly{Int64}, ::Polynomials.Poly{Int64})

julia> A
2×2 Array{Polynomials.Poly{Int64},2}:
Poly(1 + x)    Poly(2 + 2⋅x)
Poly(2 + 2⋅x)  Poly(1 + x)  

julia> A.' # Warning on no-op `transpose`
2×2 Array{Polynomials.Poly{Int64},2}:
Poly(1 + x)    Poly(2 + 2⋅x)
Poly(2 + 2⋅x)  Poly(1 + x)  

julia> A'
2×2 Array{Polynomials.Poly{Int64},2}:
Poly(1)  Poly(2)
Poly(2)  Poly(1)

Some things worth noting/discussing from the above excerpt:

Convenient ctranspose propogates to matrix operations as polyder,
transpose(p) = copy(p) would be needed to suppress the no-op transpose warning,
- here, transpose(p) = p is a faster (more efficient) solution, I know. However, to me, one should always return temporaries from mathematical operations as the result would be changed easily (and unintentionally) in cases such as p = Poly([1,1]); q = p.'; p[2] = 5 # q is unintentionally changed,
- right now, Polynomials has this problem when taking the derivative (of zeroth order), too. Yet, this is another discussion and maybe a taste of the library writer.
dot between Poly objects needs to be defined,
b'*A*b above should be re-considered (b' will have conjugated variables, whereas A*b is unconjugated).

I hope what I have listed above would make sense for you, and I hope I had not consumed a lot of your time.

Looking forward to your opinions.

A few quick comments, not touching on all you have here.

The definition of conj(p) follows from the definition conj(p)(z) = conj(p(conj(z))). This implies that for polynomials that conj(p) = sum conj(a_i) x^i. It seems that you are looking for a different property, which may be more useful.

As for the use of ctranspose for differentiation. I would be sad to see that go, but it is a pun and it seems to not play nicely with matrices of polynomials. Maybe it should be deprecated. It would be nice if there were a postfix operator in unicode that could be used.

Well, that would be cool. For now, we can define \delta<tab> = polyder(p, 1), or? (then, exporting this symbol might clash with another one coming from another package, of course).

as for the conj, I guess I am applying the correct thing, right?

\operatorname{conj}(p(x)) & = \operatorname{conj}(\sum_{n = 0}^{N} a_{n} x^{n}) \\
& = \sum_{n = 0}^{N} \operatorname{conj}(a_{n}) \operatorname{conj}(x)^{n}

(apologies for the inconvenient LaTeX up there). At least, when I thought of conj of a polynomial p, I assumed the operation was exactly as above. However, right now, Polynomials spit out the conjugated coefficients to the user and depend on the user to remember the conjugated x. It also allows for conj(p) + p operation for a polynomial p, which is again a different thing from what we would obtain from pen and pencil. Of course, the current practice is valid for all T<:Real inputs of x::T.

Maybe I have the very extreme(?) case when x is allowed to vary in the complex domain, but I just wanted to pose this question anyways.

Sorry for the confusion or wasting the time :)

I'm no expert here, but the defn I've seen is γ ∘ p ∘ γ where γ is the conjugation operator. Yours is just γ ∘ p. The defn means that conj(p)(z) = conj(p(conj(z))) in Julia, which I think the current definition satisfies.

Exactly! I was thinking of p as a mapping on the form I mentioned above, and its conjugate implies that the input (x in this case) needs to be conjugated, too. This is also obvious in the definition you have written: conj(p)(z) = conj(p(conj(z))). This means, the user is responsible for keeping track of conj(z), and this brings a bit more complications when we mix the polynomials of the form conj(p)(z) + p(z) = conj(p(conj(z))) + p(z). Here, in the resulting polynomial, we have lost z (the coefficients of p and their conjugates are added together, and we do not know which input to use).

Again, sorry for raising this issue --- apparently this is only affecting me in my special case.

I think the confusion is what exactly we are thinking p is. I'm thinking it as the function x -> p(x), not the expression p(x). The expression would have conj(p(x)) = conj(p)(conj(x)) using the current defn. of conj(p), which is what you are expecting. What is the correct interpretation for p in this package I'm not sure.

On Tue, Dec 6, 2016 at 12:13 PM, Arda Aytekin notifications@github.com wrote:

Exactly! I was thinking of p as a mapping on the form I mentioned above, and its conjugate implies that the input (x in this case) needs to be conjugated, too. This is also obvious in the definition you have written: conj(p)(z) = conj(p(conj(z))). This means, the user is responsible for keeping track of conj(z), and this brings a bit more complications when we mix the polynomials of the form conj(p)(z) + p(z) = conj(p(conj(z))) + p(z). Here, in the resulting polynomial, we have lost z (the coefficients of p and their conjugates are added together, and we do not know which input to use).

Again, sorry for raising this issue --- apparently this is only affecting me in my special case.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/Keno/Polynomials.jl/issues/85#issuecomment-265210366, or mute the thread https://github.com/notifications/unsubscribe-auth/AAZvTCEuxYMKktVzUaw0cVuFZuHFBvgbks5rFZevgaJpZM4LDwDJ .

-- John Verzani Chair, Department of Mathematics College of Staten Island, CUNY verzani@math.csi.cuny.edu

JuliaMath / Polynomials.jl