[SciPy-Dev] linalg.toeplitz and ticker #667
Warren Weckesser
warren.weckesser at enthought.com
Sat Mar 20 21:29:33 EDT 2010
josef.pktd at gmail.com wrote:
> On Mon, Mar 8, 2010 at 10:25 PM, Warren Weckesser
> <warren.weckesser at 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 at gmail.com wrote:
>>>
>>>
>>>> On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser
>>>> <warren.weckesser at 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 at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>>
>>>
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>
>>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
More information about the SciPy-Dev
mailing list