[Numpy-discussion] poly class question

Mon Oct 5 10:52:01 EDT 2009

On Mon, Oct 5, 2009 at 5:37 AM, Sebastian Walter
<sebastian.walter at gmail.com> wrote:
> On Fri, Oct 2, 2009 at 10:40 PM,  <josef.pktd at gmail.com> wrote:
>> On Fri, Oct 2, 2009 at 3:38 PM, Charles R Harris
>> <charlesr.harris at gmail.com> wrote:
>>>
>>>
>>> On Fri, Oct 2, 2009 at 12:30 PM, <josef.pktd at gmail.com> wrote:
>>>>
>>>> On Fri, Oct 2, 2009 at 2:09 PM, Charles R Harris
>>>> <charlesr.harris at gmail.com> wrote:
>>>> >
>>>> >
>>>> > On Fri, Oct 2, 2009 at 11:35 AM, Charles R Harris
>>>> > <charlesr.harris at gmail.com> wrote:
>>>> >>
>>>> >>
>>>> >> On Fri, Oct 2, 2009 at 11:33 AM, Charles R Harris
>>>> >> <charlesr.harris at gmail.com> wrote:
>>>> >>>
>>>> >>>
>>>> >>> On Fri, Oct 2, 2009 at 11:30 AM, Charles R Harris
>>>> >>> <charlesr.harris at gmail.com> wrote:
>>>> >>>>
>>>> >>>>
>>>> >>>> On Fri, Oct 2, 2009 at 11:08 AM, <josef.pktd at gmail.com> wrote:
>>>> >>>>>
>>>> >>>>> Is there a way in numpy (or scipy) to get an infinite expansion for
>>>> >>>>> the inverse of a polynomial (for a finite number of terms)
>>>> >>>>>
>>>> >>>>> np.poly1d([ -0.8, 1])**(-1)
>>>> >>>>>
>>>> >>>>> application for example the MA representation of an AR(1)
>>>> >>>>>
>>>> >>>>
>>>> >>>> Hmm, I've been working on a chebyshev class and division of a scalar
>>>> >>>> by
>>>> >>>> a chebyshev series is
>>>> >>>> expressly forbidden, but it could be included if a good interface is
>>>> >>>> proposed. Same would go for polynomials.
>>>> >>>
>>>> >>> In fact is isn't hard to get, for poly1d you should be able to
>>>> >>> multiply
>>>> >>> the series by a power of x to shift it left, then divide.
>>>> >>>
>>>> >>
>>>> >> That is, divide a power of x by the polynomial.
>>>> >>
>>>> >
>>>> > You will also need to reverse the denominator coefficients...Chuck
>>>>
>>>> That's the hint I needed. However the polynomial coefficients are then
>>>> reversed and not consistent with other polynomial operations, aren't
>>>> they?
>>>>
>>>> >>> from scipy.signal import lfilter
>>>>
>>>> >>> (np.poly1d([1, 0])**10)/np.poly1d([1, -0.8])
>>>> (poly1d([ 1.        ,  0.8       ,  0.64      ,  0.512     ,  0.4096    ,
>>>>        0.32768   ,  0.262144  ,  0.2097152 ,  0.16777216,
>>>> 0.13421773]), poly1d([ 0.10737418]))
>>>>
>>>> >>> lfilter([1], [1,-0.8], [1] + [0]*9)
>>>> array([ 1.        ,  0.8       ,  0.64      ,  0.512     ,  0.4096    ,
>>>>        0.32768   ,  0.262144  ,  0.2097152 ,  0.16777216,  0.13421773])
>>>>
>>>> >>> (np.poly1d([1, 0])**10)/np.poly1d([1, -0.8, 0.2])
>>>> (poly1d([ 1.        ,  0.8       ,  0.44      ,  0.192     ,  0.0656    ,
>>>>        0.01408   , -0.001856  , -0.0043008 , -0.00306944]),
>>>> poly1d([-0.00159539,  0.00061389]))
>>>> >>> lfilter([1], [1,-0.8, 0.2], [1] + [0]*9)
>>>> array([ 1.        ,  0.8       ,  0.44      ,  0.192     ,  0.0656    ,
>>>>        0.01408   , -0.001856  , -0.0043008 , -0.00306944, -0.00159539])
>>>>
>>>>
>>>> What I meant initally doesn't necessarily mean division of a scalar.
>>>>
>>>> >>> np.poly1d([1])/np.poly1d([-0.8, 1])
>>>> (poly1d([ 0.]), poly1d([ 1.]))
>>>>
>>>> I didn't find any polynomial division that does the expansion of the
>>>> remainder. The same problem, I think is inherited, by the
>>>> scipy.signal.lti, and it took me a while to find the usefulness of
>>>> lfilter in this case.
>>>>
>>>> If it were possible to extend the methods for the polynomial class to
>>>> do a longer expansions, it would make them more useful for arma and
>>>> lti.
>>>>
>>>> (in some areas, I'm still trying to figure out whether some
>>>> functionality is just hidden to me, or actually a limitation of the
>>>> implementation or a missing feature.)
>>>>
>>>
>>> Could you describe the sort of problems you want to solve? There are lots of
>>> curious things out there we could maybe work with. Covariances, for
>>> instance, are closely related to Chebyshev series.
>>
>> I am working on a discrete time arma process of the form
>>
>> a(L) x_t = b(L) u_t,  where L is the lag operator L^k x_t = x_(t-k)
>>
>> what I just programmed using lfilter is
>> x_t = b(L)/a(L) u_t  where  b(L)/a(L) is the impulse response function
>> or moving average representation
>>
>>  a(L)/b(L)  is the autoregressive representation
>>
>> the extension
>> a(L)(1-L)^d x_t = b(L) u_t,  where d = 0,1,2,...  (standard) or  also
>> continuous d <1 (fractional integration)
>>
>>  a(L)/b(L),  b(L)/a(L)  (1-L)^(-d)  or (1-L)^d (0<d<1)  are infinite
>> dimensional lag polynomials in the general case.
>> Initially I was looking for an easy way to do these calculation as polynomials.
>> (The fractional case (1-L)^d (0<d<1) might be a pretty special case,
>> and I just looked it up today, but is a popular model class in
>> econometrics, fractionally integrated arma processes)
>>
>> multiplication works well np.poly1d([-1, 1])*np.poly1d([-0.8, 1])
>> (with reversed poly coefficient scipy.signal I think)
>>
>> the functions in scipy.signal for lti are only for continuous time
>> processes and use poly1d under the hood, which means for example
>>
>>>>> from scipy import signal
>>>>> signal.impulse(([1, -0.8],[1]), N=10)
>>    raise ValueError, "Improper transfer function."
>> ValueError: Improper transfer function.
>>
>> while this works
>>>>> signal.impulse(([1],[1, -0.8]), N=10)
>>
>> (It's been a while since I looked inside scipy.signal.lti)
>>
>> A separate issue would be the multivariate version VARMA, or MIMO in
>> system modeling.  a(L), b(L) are matrix polynomials and x_t, u_t are
>> 1d arrays evolving in time.
>> But that is a different discussion.
>

> I'm working on something that requires truncated
> univariate/multivariate operations
> on scalar, vector and matrix polynomials. I'm pretty sure I'm doing
> something different
> than what is used in MIMO.  Still, do you have a good reference for
> implementation/complexity/stability of operations on multivariate
> (matrix) polynomials?
> Want to make sure I'm not reinventing the wheel ;)

Sorry, but I'm no help here. I was hoping to benefit from the knowledge
of others. I have used univariate and multivariate polynomials for special
cases for function approximation and time series analysis, but I know
little about the general numerical issues in this. For MIMO, I only looked
at the matlab systems toolbox, which would be an extension of scipy.signal.lti.

Since these operations are in my case usually inside an optimization loop
in a (supposedly) well behaved problem, I usually cared more about
speed and approximation errors in low order polynomials than numerical
stability.

If you invent the wheel, I would be very glad to use it.

Josef

>
>
>
>>
>> I'm not very familiar with Chebychev polynomials, the last time I
>> wanted to use them I didn't see anything about their use as a base for
>> functions in several variables and gave up. I've seen papers that use
>> them as base for functions in one variable, but I'm not doing anything
>> like this right now.
>>
>> Thanks,
>>
>> Josef
>>
>>>
>>> Chuck
>>>
>>>
>>> _______________________________________________
>>> NumPy-Discussion mailing list
>>> NumPy-Discussion at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>>
>>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>