Opened 15 years ago

Closed 12 years ago

# implement condition number for matrices

Reported by: Owned by: schilly was major sage-duplicate/invalid/wontfix linear algebra matrix, condition number Simon Spicer N/A

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

### comment:1 Changed 15 years ago by schilly

Milestone: → sage-2.10.4

### comment:2 Changed 15 years ago by dfdeshom

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:

• Computing the condition number for norms 1,2,inf. This can be done reasonably efficiently in gsl and lapack
• Computing condition estimators for other norms or when matrices get very large.

Is there any case that we would want compute a condition for norms other than 1,2 and inf? I'm curious

### Changed 12 years ago by spice

A one-line implementation of the condition number, which computes A.norm()*A.inverse().norm().

### comment:3 Changed 12 years ago by spice

Authors: → Simon Spicer matrix condition number added → N/A new → 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 12 years ago by schilly

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 12 years ago by rbeezer

Hi Simon,

Could you compare timings against #10837? (I didn't know #2512 existed.)

Rob

### comment:6 follow-up:  7 Changed 12 years ago by spice

Resolution: → duplicate needs_review → 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 12 years ago by rbeezer

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 12 years ago by mvngu

Milestone: sage-4.7 → sage-duplicate/invalid/wontfix
Note: See TracTickets for help on using tickets.