Opened 13 years ago

Closed 10 years ago

Last modified 10 years ago

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

Status badges

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)

trac_2512_matrix_condition_number.patch (2.5 KB) - added by spice 10 years ago.
A one-line implementation of the condition number, which computes A.norm()*A.inverse().norm().

Download all attachments as: .zip

Change History (9)

comment:1 Changed 13 years ago by schilly

  • Milestone set to sage-2.10.4

comment:2 Changed 13 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 10 years ago by spice

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

comment:3 Changed 10 years ago by spice

  • Authors set to Simon Spicer
  • 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 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 10 years ago by rbeezer

Hi Simon,

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

Rob

comment:6 follow-up: Changed 10 years ago by spice

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

  • Milestone changed from sage-4.7 to sage-duplicate/invalid/wontfix
Note: See TracTickets for help on using tickets.