Imports BluebitMatrix31

Module Module1
' The following code will demonstrate how to find
' the coeficients of a 3rd order polynomial 
' in order to have a least squares fit solution
' for a set of (x, y) data pairs
'
' 3rd order polynomial 
' a*1 + b*x + c*x^2 + d*x^3 = y
'
' Data set
' - x - - y -
' x0 = 2 y0 = 2
' x1 = 4 y1 = 3
' x2 = 5 y2 = 5
' x3 = 6 y3 = 7
' x4 = 8 y4 = 5
' x5 = 7 y5 = 6
'
' We need a least squares fit solution for the 
' system of equations :
'
' X * C = Y
'
' Matrix X contains the powers of xn of the dataset
' x0^0, x0^1, x0^2, x0^3
' x1^0, x1^1, x1^2, x1^3
' x2^0, x2^1, x2^2, x2^3
' x3^0, x3^1, x3^2, x3^3
' x4^0, x4^1, x4^2, x4^3
' x5^0, x5^1, x5^2, x5^3
' 
' Matrix Y contains the right hand tems y0 to yn
' 
' Matrix C contains the coeficients
'
' Using Singular Value decomposition we calculate the
' pseudoinverse of X, Xinv and then C becomes:
' 
' C = Xinv * Y

Sub Main()
Dim X As Matrix, Y As Matrix, C As Matrix

X = New MatrixClass
X.Size(6, 4)

'Entering data in the X matrix
X(0, 0) = 1 : X(0, 1) = 2 : X(0, 2) = 2 * 2 : X(0, 3) = 2 * 2 * 2
X(1, 0) = 1 : X(1, 1) = 4 : X(1, 2) = 4 * 4 : X(1, 3) = 4 * 4 * 4
X(2, 0) = 1 : X(2, 1) = 5 : X(2, 2) = 5 * 5 : X(2, 3) = 5 * 5 * 5
X(3, 0) = 1 : X(3, 1) = 6 : X(3, 2) = 6 * 6 : X(3, 3) = 6 * 6 * 6
X(4, 0) = 1 : X(4, 1) = 8 : X(4, 2) = 8 * 8 : X(4, 3) = 8 * 8 * 8
X(5, 0) = 1 : X(5, 1) = 7 : X(5, 2) = 7 * 7 : X(5, 3) = 7 * 7 * 7

Y = New MatrixClass
Y.Size(6, 1)

'Entering data for the Y matrix
Y(0, 0) = 2
Y(1, 0) = 3
Y(2, 0) = 5
Y(3, 0) = 7
Y(4, 0) = 5
Y(5, 0) = 6

'Calculating the pseudoinverse of X using SVD
Dim U As Matrix, S As Matrix, V As Matrix, Xinv As Matrix, i As Integer

U = New MatrixClass
S = New MatrixClass
V = New MatrixClass

X.SVD(U, S, V)

' Xinv = V * SInv * UInv

' Matrix S is a diagonal matrix. In order to find its
' inverse is faster to replace its diagonal elements by
' their reciprocals
For i = 0 To 3
S(i, i) = 1 / S(i, i)
Next


' Matrix U is an orthogonal matrix - its inverse equals
' its transpose

' Finally the X pseudoinverse becomes :
Xinv = V.Times(S).Times(U.Transpose)

' and the C matrix containg the coeficients :

C = Xinv.Times(Y)
'now C(0,0) contains a
' C(1,0) contains b
' C(2,0) contains c
' C(3,0) contains d

' In order to generate y based on C coeficients
' we just need to multiply matrix X by C
Dim nY As Matrix = X.Times(C)

' calculating DY = Y - nY
Dim dY As Matrix = Y.Minus(nY)

' calculating the sum of squares
Dim sumDY2 As Double = dY.ColsDotProduct(0, 0)

'Writing down the results:
Console.WriteLine("Dataset:")
Console.WriteLine("- x - - y -")
For i = 0 To 5
Console.WriteLine(String.Format(" x{0}={1} y{0}={2}", i, X(i, 1), Y(i, 0)))
Next

Console.WriteLine()
Console.WriteLine("3rd order polynomial coeficients.")
Console.WriteLine(String.Format("a= {0:+0.0000;-0.0000}", C(0, 0)))
Console.WriteLine(String.Format("a= {0:+0.0000;-0.0000}", C(1, 0)))
Console.WriteLine(String.Format("a= {0:+0.0000;-0.0000}", C(2, 0)))
Console.WriteLine(String.Format("a= {0:+0.0000;-0.0000}", C(3, 0)))

Console.WriteLine()
Console.WriteLine("3rd order polynomial.")
Console.WriteLine(String.Format("{0: + 0.0000}{1: + 0.0000; - 0.0000}*x{2: + 0.0000; - 0.0000}*x^2{3: + 0.0000; - 0.0000}*x^3", C(0, 0), C(1, 0), C(2, 0), C(3, 0)))

Console.WriteLine()
Console.WriteLine(" = Y values =")
Console.WriteLine("Observed Estimated Error")
For i = 0 To 5
Console.WriteLine(String.Format(" {0:+0.0000;-0.0000} {1:+0.0000;-0.0000} {2:+0.0000;-0.0000}", Y(i, 0), nY(i, 0), dY(i, 0)))
Next
Console.WriteLine()
Console.WriteLine(String.Format("Sum of error squares : {0:+0.0000;-0.0000}", sumDY2))

Console.Read()
End Sub

End Module