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
>> 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>