<br><br>On Tuesday, June 26, 2012, Charles R Harris  wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br><br><div class="gmail_quote">On Tue, Jun 26, 2012 at 3:42 PM, Matthew Brett <span dir="ltr"><<a>matthew.brett@gmail.com</a>></span> wrote:<br>
<blockquote style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi,<br>
<div><div><br>
On Mon, Jun 18, 2012 at 3:50 PM, Matthew Brett <<a>matthew.brett@gmail.com</a>> wrote:<br>
> Hi,<br>
><br>
> On Sun, Jun 17, 2012 at 7:22 PM, Charles R Harris<br>
> <<a>charlesr.harris@gmail.com</a>> wrote:<br>
>><br>
>><br>
>> On Sat, Jun 16, 2012 at 2:33 PM, Matthew Brett <<a>matthew.brett@gmail.com</a>><br>
>> wrote:<br>
>>><br>
>>> Hi,<br>
>>><br>
>>> On Sat, Jun 16, 2012 at 8:03 PM, Matthew Brett <<a>matthew.brett@gmail.com</a>><br>
>>> wrote:<br>
>>> > Hi,<br>
>>> ><br>
>>> > On Sat, Jun 16, 2012 at 10:40 AM, Nathaniel Smith <<a>njs@pobox.com</a>> wrote:<br>
>>> >> On Fri, Jun 15, 2012 at 4:10 AM, Charles R Harris<br>
>>> >> <<a>charlesr.harris@gmail.com</a>> wrote:<br>
>>> >>><br>
>>> >>><br>
>>> >>> On Thu, Jun 14, 2012 at 8:06 PM, Matthew Brett<br>
>>> >>> <<a>matthew.brett@gmail.com</a>><br>
>>> >>> wrote:<br>
>>> >>>><br>
>>> >>>> Hi,<br>
>>> >>>><br>
>>> >>>> I noticed that numpy.linalg.matrix_rank sometimes gives full rank for<br>
>>> >>>> matrices that are numerically rank deficient:<br>
>>> >>>><br>
>>> >>>> If I repeatedly make random matrices, then set the first column to be<br>
>>> >>>> equal to the sum of the second and third columns:<br>
>>> >>>><br>
>>> >>>> def make_deficient():<br>
>>> >>>>    X = np.random.normal(size=(40, 10))<br>
>>> >>>>    deficient_X = X.copy()<br>
>>> >>>>    deficient_X[:, 0] = deficient_X[:, 1] + deficient_X[:, 2]<br>
>>> >>>>    return deficient_X<br>
>>> >>>><br>
>>> >>>> then the current numpy.linalg.matrix_rank algorithm returns full rank<br>
>>> >>>> (10) in about 8 percent of cases (see appended script).<br>
>>> >>>><br>
>>> >>>> I think this is a tolerance problem.  The ``matrix_rank`` algorithm<br>
>>> >>>> does this by default:<br>
>>> >>>><br>
>>> >>>> S = spl.svd(M, compute_uv=False)<br>
>>> >>>> tol = S.max() * np.finfo(S.dtype).eps<br>
>>> >>>> return np.sum(S > tol)<br>
>>> >>>><br>
>>> >>>> I guess we'd we want the lowest tolerance that nearly always or<br>
>>> >>>> always<br>
>>> >>>> identifies numerically rank deficient matrices.  I suppose one way of<br>
>>> >>>> looking at whether the tolerance is in the right range is to compare<br>
>>> >>>> the calculated tolerance (``tol``) to the minimum singular value<br>
>>> >>>> (``S.min()``) because S.min() in our case should be very small and<br>
>>> >>>> indicate the rank deficiency. The mean value of tol / S.min() for the<br>
>>> >>>> current algorithm, across many iterations, is about 2.8.  We might<br>
>>> >>>> hope this value would be higher than 1, but not much higher,<br>
>>> >>>> otherwise<br>
>>> >>>> we might be rejecting too many columns.<br>
>>> >>>><br>
>>> >>>> Our current algorithm for tolerance is the same as the 2-norm of M *<br>
>>> >>>> eps.  We're citing Golub and Van Loan for this, but now I look at our<br>
>>> >>>> copy (p 261, last para) - they seem to be suggesting using u * |M|<br>
>>> >>>> where u = (p 61, section 2.4.2) eps /  2. (see [1]). I think the<br>
>>> >>>> Golub<br>
</div></div></blockquote><div><br>I'm fine with that, and agree that it is likely to lead to fewer folks wondering why Matlab and numpy are different. A good explanation in the function documentation would be useful.<br>

<br>Chuck<br></div><br></div></blockquote><div><br></div><div>One potential problem is that it implies that it will always be the same as any version of matlab's tolerance.  What if they change it in a future release? How likely are we to even notice?</div>
<div><br></div><div>Cheers!</div><div>Ben Root <span></span></div>