#2512 closed enhancement (duplicate)
implement condition number for matrices
Reported by: | schilly | Owned by: | was |
---|---|---|---|
Priority: | major | Milestone: | sage-duplicate/invalid/wontfix |
Component: | linear algebra | Keywords: | matrix, condition number |
Cc: | Merged in: | ||
Authors: | Simon Spicer | Reviewers: | |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description
implement condition number (wikipedia) for matrices.
something like:
def condition(m,p=2): return norm(m.inverse(),p) * norm(m,p)
depends on #1763
Attachments (1)
Change History (9)
comment:1 Changed 13 years ago by
- Milestone set to sage-2.10.4
comment:2 Changed 13 years ago by
Changed 10 years ago by
A one-line implementation of the condition number, which computes A.norm()*A.inverse().norm().
comment:3 Changed 10 years ago by
- Keywords matrix condition number added
- Report Upstream set to N/A
- Status changed from new to needs_review
The patch implements a very simple .condition_number() method, based on a matrix's .norm() method. As such it only works for p = 1,2,Infinity or the Frobenius norm.
Suggestions/comments welcome. Written on sage 4.6.2.
comment:4 Changed 10 years ago by
wow, this still exists. The doctest looks fine, is it possible to indent the options for the argument? Apart from that I strongly suggest to use numpy for that. They solve the 2-norm special case via svd and all the others might be faster. They also have a -infinity case.
sage: from numpy.linalg.linalg import cond sage: A = matrix([[1,2,4],[5,3,9],[7,8,6]]) sage: cond(A) 13.139632629120618 sage: cond(A, 1) 22.88636363636364 sage: cond(A, 'fro') 14.21690371493278 sage: import numpy as np sage: cond(A, np.inf) 19.090909090909093 sage: cond(A, -np.inf) 2.5454545454545454
complex also works
sage: A = matrix([[1+2j, 1+3j], [1+1j, 0.5-0.5j]]) sage: cond(A, np.inf) 4.2200687516284452 sage: cond(A) 3.2255049266776936
comment:5 Changed 10 years ago by
comment:6 follow-up: ↓ 7 Changed 10 years ago by
- Resolution set to duplicate
- Status changed from needs_review to closed
Turns out Rob Beezer has implemented condition numbers by directly wrapping the NumPy? cond() command in patch #10837.
Some timings:
Patch 2512:
sage: A = matrix([[1,2,4],[5,3,9],[7,8,6]]) sage: timeit('A.condition_number(2)') 125 loops, best of 3: 1.96 ms per loop sage: timeit('A.condition_number(Infinity)') 625 loops, best of 3: 1.31 ms per loop sage: timeit('A.condition_number("frob")') 625 loops, best of 3: 949 µs per loop sage: A = MatrixSpace(CC,50).random_element() sage: timeit('A.condition_number(2)') 5 loops, best of 3: 389 ms per loop sage: timeit('A.condition_number(1)') 5 loops, best of 3: 380 ms per loop sage: timeit('A.condition_number("frob")') 5 loops, best of 3: 358 ms per loop
Patch 10837:
sage: A = matrix(RDF,[[1,2,4],[5,3,9],[7,8,6]]) sage: timeit('A.condition(2)') 625 loops, best of 3: 282 µs per loop sage: timeit('A.condition(Infinity)') 625 loops, best of 3: 123 µs per loop sage: timeit('A.condition("frob")') 625 loops, best of 3: 203 µs per loop sage: A = MatrixSpace(CDF,50).random_element() sage: timeit('A.condition(1)') 625 loops, best of 3: 968 µs per loop sage: timeit('A.condition(2)') 125 loops, best of 3: 2.97 ms per loop sage: timeit('A.condition("frob")') 625 loops, best of 3: 942 µs per loop
I vote we go for Rob's patch.
comment:7 in reply to: ↑ 6 Changed 10 years ago by
I vote we go for Rob's patch.
Yep, looks like 3 to 10 times faster. I might steal your code for a verification doctest on #10837.
Thanks for the timing tests.
comment:8 Changed 10 years ago by
- Milestone changed from sage-4.7 to sage-duplicate/invalid/wontfix
Condition numbers for any norm can be expensive, since they involve computing the inverse of a matrix. If our matrix is singular, this is even worse. I would propose:
Is there any case that we would want compute a condition for norms other than 1,2 and inf? I'm curious