On 8/6/15, Waldek Hebisch <hebi...@math.uni.wroc.pl> wrote:
> Abhinav Baid wrote:
>>
>> On 8/5/15, Waldek Hebisch <hebi...@math.uni.wroc.pl> wrote:
>> >>
>> >> Just thought that I'd make a small post. I've not been able to work on
>> >> this much for the past week as college reopened. Waldek, have you gone
>> >> through the current code? Is it fine, or could you please point out
>> >> the errors?
>> >
>> > In general it look OK. However, I noticed that you use
>> > the same bounds for operator and its adjoint. Generalized
>> > exponents of adjoint of L are related, but usually different
>> > than generalized exponents of L. So, I would expect different
>> > bound for adjoint. Using too small bound may lead to missing
>> > a solution.
>>
>> Okay, I use a different bound for adjoint now. [1]
>
> Good. AFAICS there is still serious problem: you pass too small
> number of coefficients to Hermite-Pade solver. The rule is:
> if you want to get back n coefficients you need to have at
> least n equations. However, the number of equations is
> related to order of approximation. In scalar Hermite-Pade
> problem number of equations is the same as order of approximation.
> So to get n coefficients you need approximation of order n,
> which means that you need n terms for _each_ coefficient of
> local factor. For example, if we seek factor of order 1
> and degree bound is 10, that we need to determine 22 numbers
> (11 for each coefficient of the factor). For this we need
> 22 terms of series. Note that missing coefficients are
> effectively replaced by 0 (since we pass (list of vectors of)
> polynomials to Hermite-Pade solver). So we are solving
> wrong Hermite-Pade problem and conseqently we should
> expect wrong solutions.
>
> Additional remarks:
> - If there is singularity with exponential
> part of multiplicity 1, then failure of van Hoej algorithm
> for this exponential part applied to operator and its
> adjount means that operator is irreducible. In such
> case there is no need to look at other singularities
> or even other exponnetial parts. And if our implementation
> fails to find factor for reducible operator this is bug
> which should be fixed.
Yes, but how do I check this and appropriately call factor_minmult1?
(which is removed for now). The earlier check was using
skipped_factors but it was wrong. I think there should be a condition
in the loop after bounds are computed, but am not sure. Also, I have
fixed the code according to the other remarks.
> - In try_factorization2 you perform too many calls to
> Hermite-Pade solver. Normally a single call per order
> of factor should be enough (for possible exceptions
> see remarks at the end)
> - You have some redundant code, for example you call
> 'truncate' and then 'convertUTStoUP' which in turn
> callse 'UTS2UP'. But 'UTS2UP' ignores unneeded terms
> of the series, so call to 'truncate' is not needed.
> - Could you try to use better variable names. 'i', 'j', 'k'
> usually mean integer variables. Unless there is strong
> reason (like trying to follow notation in a paper) please
> do not use then for non-integer variables.
>
>> > Also, it is not clear if you can just pick
>> > a basis element for a Hermite-Pade problem: columns of solution
>> > matrix span solutions space. To prove that operator is
>> > irreducible you need to check that no element of solution
>> > space leads to a factor. Unfortunetely, currently 'fffg'
>> > routine will always return several solutions. More precisely,
>> > 'fffg' only accepts bounds that lead to several solutions.
>> > I will try to modify 'fffg' so that it accepts bounds leading
>> > to smaller number of solutions (hopefully a unique solution).
>> >
>>
>> So, as the solutions are spanned by the columns of the matrix, I need
>> to find coefficients, not all zero, for each of them such that the
>> corresponding linear combination gives a dot product of zero with lv.
>> Is there a relevant function for this? I looked at expressIdealMember,
>> but it returns the trivial list of zeros.
>
> I have commited a new version of fffg which accepts requested
> order of solution as an argument. Normally if requested order
> is large enoung you will get single solution. If you look
> for factor starting from factors of low order, then if you
> take sufficiently many terms (that is order of Hermite-Pade
> problem is large enough) you will get single solution. ATM
> I do not know if there is a simple way to find out which
> order is large enough. One way could be to increase number
> of terms in Hermite-Pade when there are multiple solutions
> and stop once we get only one.
>
> Below is a modified version of 'hp_solve' which accepts
> requested order of solution as an argument.
>
> Note: New version of 'generalInterpolation' computes
> defects of solutions, that is margin between degree
> of solution and degree bound. If defect is negative
> then given "solution" is really not a solution because
> it does not satisfy degree bounds. 'hp_solve' rejects
> them, so result of 'hp_solve' consists of true solutions
> (in case there are no solution return value is matrix
> with 0 columns).
>
>
> ----------------<cut here>-------------------
>
> )abbrev package VHPSOLV VectorHermitePadeSolver
> VectorHermitePadeSolver : Exports == Implementation where
> F ==> Expression(Integer)
> SUP ==> SparseUnivariatePolynomial(F)
> Exports ==> with
> hp_solve : (List Vector(SUP), List(NonNegativeInteger),
> NonNegativeInteger) -> Matrix SUP
> +++ hp_solve(lv, eta, K) solves Hermite-Pade problem with degree
> +++ bound eta up to order K.
> Implementation ==> add
>
> FFFG ==> FractionFreeFastGaussian(F, SUP)
>
> power_action(m : NonNegativeInteger
> ) : ((NonNegativeInteger, NonNegativeInteger, SUP) -> F) ==
> (k, l, g) +-> DiffAction(k, m*l, g)$FFFG
>
> hp_solve(lv, eta, K) ==
> print("hp_solve"::OutputForm)
> print(eta::OutputForm)
> m := #first(lv)
> lpp : List(SUP) := []
> for v in lv repeat
> #v ~= m =>
> error "hp_solve: vectors must be of the same length"
> pp : SUP := 0
> for i in 1..m repeat
> pp1 := multiplyExponents(v(i), m)
> pp := pp + monomial(1, (i-1)::NonNegativeInteger)$SUP*pp1
> lpp := cons(pp, lpp)
> lpp := reverse!(lpp)
> vd : Vector Integer := vector([ei::Integer for ei in eta])
> C : List(F) := [0 for i in 1..K]
> res1 := generalInterpolation(C, power_action(m), vector(lpp),
> vd, K)$FFFG
> print(vd::OutputForm)
> n := #vd
> n1 : NonNegativeInteger := 0
> for i in 1..n | vd(i) >= 0 repeat
> n1 := n1 + 1
> print(n1::OutputForm)
> res := new(n, n1, 0)$Matrix(SUP)
> i1 : NonNegativeInteger := 0
> for i in 1..n | vd(i) >= 0 repeat
> i1 := i1 + 1
> for j in 1..n repeat
> res(j, i1) := res1(j, i)
> res
>
> ------------------<cut here>-------------------
> --
> Waldek Hebisch
> hebi...@math.uni.wroc.pl
>
> --
> You received this message because you are subscribed to a topic in the
> Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/fricas-devel/3MYCZGOyIOI/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> fricas-devel+unsubscr...@googlegroups.com.
> To post to this group, send email to fricas-devel@googlegroups.com.
> Visit this group at http://groups.google.com/group/fricas-devel.
> For more options, visit https://groups.google.com/d/optout.
>
--
Abhinav.
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to fricas-devel+unsubscr...@googlegroups.com.
To post to this group, send email to fricas-devel@googlegroups.com.
Visit this group at http://groups.google.com/group/fricas-devel.
For more options, visit https://groups.google.com/d/optout.
