[WIP] Adding norm and cond function for issue #436

[arafatkatze]

I am picking up the work done from pr #400 and #204.
This Pr is to solve the issue #436 and #470 using cond function and also to add norm function as per issue #54 .

Condition Number

The condition number of a matrix measures the sensitivity of the solution of a system of linear equations to errors in the data. It gives an indication of the accuracy of the results from matrix inversion and the linear equation solution. Values of cond(X) and cond(X,p) near 1 indicate a well-conditioned matrix.

A very high value of condition number indicates instability.

Example for singular matrix

[1] pry(main)> a = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0], dtype: :int32)
=> 
[
  [ 1,  2, 1]   [-2, -3, 1]   [ 3,  5, 0] ]
[2] pry(main)> a.det
=> 0
[3] pry(main)> a.invert
=> 
[
  [ 7.505999378950824e+15, -7.505999378950824e+15, -7.505999378950824e+15]
  [-4.503599627370494e+15,  4.503599627370494e+15,     4503599627370494.5]
  [    1501199875790165.2,    -1501199875790164.2,    -1501199875790164.8]
]
[4] pry(main)> a.cond
=> 1.3510798882111485e+17

Note that the very high cond value is good indicator of instability. 

Example for non-singular matrix

Now taking another non-singular matrix 
[10] pry(main)> a = NMatrix.new([3, 3], [909325235,  201242,  1213, -20213, -339921,  1421,  43,      55, 99123], dtype: :int32)
=> 
[
  [909325235,  201242,  1213]   [   -20213, -339921,  1421]   [       43,      55, 99123] ]
[11] pry(main)> a.det
=> -30638462559409565696
[12] pry(main)> a.invert
=> 
[
  [   1.09973101204623e-09,     6.5106543816683e-10, -2.279125637084372e-11]
  [ -6.539604583992586e-11, -2.9418919125582587e-06,  4.217495165752539e-08]
  [-4.4078216959527015e-13,  1.6320738817112382e-09, 1.0088452542373443e-05]
]
[13] pry(main)> a.cond
=> 9212.261157823907

Note here that this cond value is much less than 1.3510798882111485e+17 as in previous example
Note: The exact value of condition number varies with implementation as mentioned in another comment below.

TODO :

• Add more tests for norm function with support for all data types
• Check for overflow errors
• Checking if the results for non-singular matrices have cond function value much less than a threshold
• Add tests for cond function
• Add warning for ill-conditioned matrices

This pr is not complete and there are mistakes(as mentioned in a comment below) and I will work on them. I can understand if I might be missing something,and it would be very good if someone could guide me.

[arafatkatze]

I am experiencing this issue that inverse is implemented differently when using plugins

These are a few examples on a same matrix under different plugins.

[1] pry(main)> b = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0], dtype: :float64)
=> 
[
  [ 1.0,  2.0, 1.0]   [-2.0, -3.0, 1.0]   [ 3.0,  5.0, 0.0] ]
[2] pry(main)> b.invert
=> 
[
  [ 7.505999378950824e+15, -7.505999378950824e+15, -7.505999378950824e+15]
  [-4.503599627370494e+15,  4.503599627370494e+15,     4503599627370494.5]
  [    1501199875790165.2,    -1501199875790164.2,    -1501199875790164.8]
]
[3] pry(main)> b.cond
=> 1.3510798882111485e+17

lapacke

[1] pry(main)> require 'nmatrix/lapacke'
=> true
[2] pry(main)> b = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0], dtype: :float64)
=> 
[
  [ 1.0,  2.0, 1.0]   [-2.0, -3.0, 1.0]   [ 3.0,  5.0, 0.0] ]
[3] pry(main)> b.invert
=> 
[
  [-3.752999689475412e+15,  3.752999689475412e+15, 3.752999689475412e+15]
  [ 2.251799813685247e+15, -2.251799813685247e+15,   -2251799813685246.8]
  [    -750599937895081.6,      750599937895082.6,     750599937895082.2]
]
[4] pry(main)> b.cond
=> 6.755399441055742e+16

Atlas

[1] pry(main)> require 'nmatrix/atlas'
=> true
[2] pry(main)> b = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0], dtype: :float64)
=> 
[
  [ 1.0,  2.0, 1.0]   [-2.0, -3.0, 1.0]   [ 3.0,  5.0, 0.0] ]
[3] pry(main)> b.invert
=> 
[
  [ 2.0,  0.5,  0.5]   [-3.0,  2.5, -0.2]   [ 5.0, -2.5, -0.0] ]
[4] pry(main)> b.cond
=> 100.0

The examples given above indicate that on using ATLAS for calculating the inverse and the condition number am getting incorrect results as the condition number is very low.

Note that the exact value for condition number is not important the only requirement is
that it must be very high to indicate ill conditioned matrix.

The value of condition number is different under different implementations for example

Numpy

from numpy import *
import numpy as np
a = np.matrix([   [ 1.0,  2.0, 1.0] ,  [-2.0, -3.0, 1.0] ,  [ 3.0,  5.0, 0.0]  ])
print np.linalg.det(a)  
#1.33226762955e-15  
print np.linalg.cond(a,1)
#6.75539944106e+16

Sympy

from sympy import *
a = Matrix(3,3,[ 1.0,  2.0, 1.0,-2.0, -3.0, 1.0, 3.0,  5.0, 0.0] )
print a.det()
# 0
b= a.inv()
# Raises error that the matrix is non invertible hence norm cannot be calculated

Octave

octave:24> a =  [1.0,2.0,1.0;-2.0,-3.0,1.0;3.0,5.0,0]
a =

   1   2   1
  -2  -3   1
   3   5   0

octave:25> det(a)
ans =    8.3267e-16
octave:26> cond(a,1)
ans =    1.0809e+17

[wlevine]

Can you file a bug on the atlas implementation? I'll try to figure out what's going on.

[translunar]

Just as a quick note, we probably want an rcond function in addition to cond to prevent overflow and for consistency with other libraries.