[SciPy-Dev] linalg.toeplitz and ticker #667

Warren Weckesser warren.weckesser@enthought....
Sat Mar 20 20:29:33 CDT 2010


josef.pktd@gmail.com wrote:
> On Mon, Mar 8, 2010 at 10:25 PM, Warren Weckesser
> <warren.weckesser@enthought.com> wrote:
>   
>> Two more questions about the toeplitz API:
>>     
>
> I did a (partial) search in trac and on the mailing lists and didn't
> find any information
>
>   


I have submitted a patch, attached to ticket #667.   The patch changes
the behavior of toeplitz() and hankel(), and adds a new function, 
circulant().
The patch also adds several tests of these functions.  I can commit the
patch, but I'd like to know if y'all agree with the changes first.

Warren


>   
>> The signature is toeplitz(c, r=None).  In the current implementation,
>> if either c or r is a scalar, then c is returned.  This results in the
>> following:
>>
>> In [1]: from scipy.linalg import toeplitz
>>
>> In [2]: toeplitz(1, [1,2,3,4])
>> Out[2]: 1
>>
>> In [3]: toeplitz([1], [1,2,3,4])
>> Out[3]: array([[1, 2, 3, 4]])
>>
>> I would expect these to produce the same result--the second result,
>> in fact.
>>
>> Similarly, the following is surprising:
>>
>> In [4]: toeplitz([1,2,3,4], 1)
>> Out[4]: [1, 2, 3, 4]
>>
>> In [5]: toeplitz([1,2,3,4], [1])
>> Out[5]:
>> array([[1],
>>       [2],
>>       [3],
>>       [4]])
>>
>> I think it makes more sense for toeplitz to simply treat a scalar
>> as a 1D array of length 1.  Moreover, it seems to me that toeplitz
>> should *always* return a 2D array.
>>     
>
> numpy/scipy is not matlab, so the policy is still a bit vague
>
>   
>>>> np.kron(5,3)
>>>>         
> 15
>   
>>>> np.kron([5],[3])
>>>>         
> array([15])
>   
>>>> np.kron([5],[3]).shape
>>>>         
> (1,)
>   
>>>> scipy.linalg.hankel(1)
>>>>         
> 1
>   
>>>> scipy.linalg.hankel(1,2)
>>>>         
> 1
>   
>>>> np.triu(3)
>>>>         
> Traceback (most recent call last):
>   File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py",
> line 442, in triu
>     out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m)
> IndexError: tuple index out of range
>   
>>>> np.triu([3])
>>>>         
> Traceback (most recent call last):
>   File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py",
> line 442, in triu
>     out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m)
> IndexError: tuple index out of range
>   
>>>> np.triu([[3]])
>>>>         
> array([[3]])
>
>   
>>>> scipy.linalg.kron([-1j], [1.j])
>>>>         
> Traceback (most recent call last):
>   File "\Programs\Python25\Lib\site-packages\scipy\linalg\basic.py",
> line 835, in kron
> AttributeError: 'list' object has no attribute 'flags'
>
> and I filed a ticket for scipy.linalg.block_diag because there I can
> see a use for scalar to 2d.
>
> I don't have a strong opinion about always 2d result, but a bit more
> consistency in numpy/scipy would be useful.
>
>
>   
>> My second question: does anyone know why toeplitz uses
>> asarray_chkfinite()?  If, in general, numpy arrays can hold the
>> values nan and inf, why should this function care if its values
>> are finite?  Any objections to not using asarray_chkfinite()?
>>     
>
> I agree, I don't think it's the job of toeplitz or hankle to throw an
> exception if I want to have a inf or nan in my matrix.
>
>   
>> Warren
>>
>>
>>
>> Warren Weckesser wrote:
>>     
>>> josef.pktd@gmail.com wrote:
>>>
>>>       
>>>> On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser
>>>> <warren.weckesser@enthought.com> wrote:
>>>>
>>>>
>>>>         
>>>>> Warren Weckesser wrote:
>>>>>
>>>>>
>>>>>           
>>>>>> I am going to update linalg.toeplitz to fix the problem reported in
>>>>>> ticket #667, but I have some questions about the desired behavior of
>>>>>> the function when either `c` (the first column) or `r` (the first row)
>>>>>> is complex.  The docstring does not say anything about complex values,
>>>>>> but the code does some undocumented and unexpected things in this
>>>>>> case, one example of which is the bug reported in the ticket #667.
>>>>>>
>>>>>> I would like to fix this, but I would first like opinions on exactly
>>>>>> how it *should* handle complex arguments.  I *think* the idea is that
>>>>>> if `c` is complex and `r` is not given, then a Hermitian Toeplitz
>>>>>> matrix is created.  Is that correct?  Note that this requires that
>>>>>> c[0] be real, since c[0] is the value that will be in the diagonal of
>>>>>> the matrix.
>>>>>>
>>>>>>
>>>>>>
>>>>>>             
>>> It looks like the behavior was copied from Matlab:
>>>
>>> http://www.mathworks.com/access/helpdesk/help/techdoc/ref/toeplitz.html
>>>
>>>
>>>       
>>>>>> While I am fixing this, I would like to remove the printed warning
>>>>>> that occurs when r[0] != c[0].  I would like to change this to raise
>>>>>> a ValueError--that is, adopt the philosophy that either the arguments
>>>>>> are correct (and the code handles them as [not yet, but soon to be]
>>>>>> documented in the docstring, with no warnings needed), or the
>>>>>> arguments are wrong and an exception should be raised.
>>>>>>
>>>>>> This means that when `r` is not given, an exception would be raised
>>>>>> if c[0] is not real.
>>>>>>
>>>>>> A very different way to handle the case of `r` not being given is to
>>>>>> assume it means a square matrix should be created with r[1:] = 0
>>>>>> (i.e. zeros in all the upper off-diagonals).  This avoids the whole
>>>>>> issue of complex numbers and conjugates, but it is also a more
>>>>>> drastic change.
>>>>>>
>>>>>>
>>>>>>             
>>>> a lot of code relies on toeplitz of a single argument to build the
>>>> symmetric matrix.
>>>>
>>>>
>>>>
>>>>         
>>> OK, thanks--that's good to know.
>>>
>>>
>>>       
>>>>>>             
>>>>> By the way, that is how linalg.hankel handles the case where `r` is not
>>>>> given.
>>>>>
>>>>>
>>>>>           
>>>> ??
>>>>
>>>>
>>>>
>>>>         
>>> I'm not sure what your question marks mean.  What I meant by "that is
>>> how linalg.hankel handles [it]" was that when `r` is not given, hankel
>>> uses zeros to fill in the bottom row.  Your examples demonstrate this.
>>> My point is moot, since if there is "a lot of code" that relies on
>>> creating a symmetric/hermitian matrix, it would be a bad idea to make
>>> such a significant change to the API.
>>>
>>>
>>>       
>>>>>>> import scipy.linalg
>>>>>>> scipy.linalg.toeplitz(np.arange(4))
>>>>>>>
>>>>>>>
>>>>>>>               
>>>> array([[0, 1, 2, 3],
>>>>        [1, 0, 1, 2],
>>>>        [2, 1, 0, 1],
>>>>        [3, 2, 1, 0]])
>>>>
>>>>
>>>>         
>>>>>>> scipy.linalg.hankel(np.arange(4))
>>>>>>>
>>>>>>>
>>>>>>>               
>>>> array([[ 0.,  1.,  2.,  3.],
>>>>        [ 1.,  2.,  3.,  0.],
>>>>        [ 2.,  3.,  0.,  0.],
>>>>        [ 3.,  0.,  0.,  0.]])
>>>>
>>>>
>>>>         
>>>>>>> scipy.linalg.hankel(np.arange(4)+1.j)
>>>>>>>
>>>>>>>
>>>>>>>               
>>>> array([[ 0.+1.j,  1.+1.j,  2.+1.j,  3.+1.j],
>>>>        [ 1.+1.j,  2.+1.j,  3.+1.j,  0.+0.j],
>>>>        [ 2.+1.j,  3.+1.j,  0.+0.j,  0.+0.j],
>>>>        [ 3.+1.j,  0.+0.j,  0.+0.j,  0.+0.j]])
>>>>
>>>>
>>>>         
>>>>>>> scipy.linalg.toeplitz(np.arange(4)+1.j)
>>>>>>>
>>>>>>>
>>>>>>>               
>>>> Warning: column and row values don't agree; column value used.
>>>> array([[ 0.+1.j,  1.+1.j,  2.+1.j,  3.+1.j],
>>>>        [ 1.-1.j,  0.+1.j,  1.+1.j,  2.+1.j],
>>>>        [ 2.-1.j,  1.-1.j,  0.+1.j,  1.+1.j],
>>>>        [ 3.-1.j,  2.-1.j,  1.-1.j,  0.+1.j]])
>>>>
>>>>
>>>> requiring the diagonal to be real might be too tough, I'm not using
>>>> complex toeplitz, but for other cases, I often have "almost real"
>>>> values, up to numerical imprecision.
>>>>
>>>>
>>>>
>>>>         
>>> This could be dealt with like Matlab does: don't worry about it, and
>>> just use c[1:] to fill in r[1:].  If c[0] is has nonzero imaginary part,
>>> then the matrix is not Hermitian--but that's OK, since that is what the
>>> user gave to the function.
>>>
>>> Your example illustrates one of the bugs: note that, except for the
>>> first element, the first column is the *conjugate* of the first
>>> argument.  The docstring says that `c` is the first column (and the
>>> optional second argument is the first row).  So it is the elements in
>>> the first row starting in the second column that should be conjugated.
>>> It also prints that annoying warning.
>>>       
>
> That's how I interpreted the bug
>
>
>   
>>>       
>>>>>>> scipy.linalg.toeplitz(np.arange(4)*1.j,-np.arange(4)*1.j)
>>>>>>>
>>>>>>>
>>>>>>>               
>>>> array([[ 0.+0.j,  0.-1.j,  0.-2.j,  0.-3.j],
>>>>        [ 0.+1.j,  0.+0.j,  0.-1.j,  0.-2.j],
>>>>        [ 0.+2.j,  0.+1.j,  0.+0.j,  0.-1.j],
>>>>        [ 0.+3.j,  0.+2.j,  0.+1.j,  0.+0.j]])
>>>>
>>>> also for real values  r[0] != c[0] might not be satisfied if there is
>>>> numerical imprecision and they are only almost equal.
>>>>
>>>>
>>>>
>>>>         
>>> True.  So if the requirement that r[0] == c[0] is strict, a user who
>>> knows the first values are "close enough", but maybe not exactly equal,
>>> would have to do something like:
>>>
>>> r[0] = c[0]
>>> t = toeplitz(c, r)
>>>
>>> I could live with that, but not everyone will agree.
>>>
>>> Instead, the function could require that the values be "close", using
>>> something like numpy.allclose(), but how close is close enough?  Should
>>> the function then have `atol` and `rtol` arguments to use for the
>>> comparison?  This feels like overkill, and is not my preference.
>>>
>>> So that leaves the behavior used by Matlab: use c[0], and ignore r[0].
>>> But I don't like the printed warning (or any logged warning) when c[0]
>>> != r[0].  If the program is behaving correctly, and as documented, no
>>> warning should be necessary.
>>>       
>
> raising a proper Warning might be useful and avoid the almost equal issue
>
> I think, all the uses I have seen use only a single argument to build
> the symmetric toeplitz matrix.
>
> Josef
>
>
>   
>>> (If I could start from scratch, I would simply define the `r` argument
>>> to give the values from the second column onward, and so completely
>>> avoid the issue.)
>>>
>>> Warren
>>>
>>> _______________________________________________
>>> SciPy-Dev mailing list
>>> SciPy-Dev@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>>
>>>       
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>
>>     
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>   



More information about the SciPy-Dev mailing list