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

josef.pktd@gmai... josef.pktd@gmai...
Sat Mar 20 23:34:43 CDT 2010


On Sat, Mar 20, 2010 at 9:29 PM, Warren Weckesser
<warren.weckesser@enthought.com> wrote:
> 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.

fine with me

Josef

>
> 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
>>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>


More information about the SciPy-Dev mailing list