Structured Linear Equations in F# QuickStart Sample
Illustrates how to solve systems of simultaneous linear equations that have special structure in F#.
View this sample in: C# Visual Basic IronPython
// Illustrates solving symmetrical and triangular systems
// of simultaneous linear equations using classes
// in the Numerics.NET.LinearAlgebra namespace of Numerics.NET.
#light
open System
open Numerics.NET
// The structured matrix classes reside in the
// Numerics.NET.LinearAlgebra namespace.
open Numerics.NET.LinearAlgebra
// The license is verified at runtime. We're using a 30 day trial key here.
// For more information, see:
// https://numerics.net/trial-key
let licensed = Numerics.NET.License.Verify("64542-18980-57619-62268")
// To learn more about solving general systems of
// simultaneous linear equations, see the
// LinearEquations QuickStart Sample.
//
// The methods and classes available for solving
// structured systems of equations are similar
// to those for general equations.
//
// Triangular systems and matrices
//
printfn "Triangular matrices:"
// For the basics of working with triangular
// matrices, see the TriangularMatrices QuickStart
// Sample.
//
// Let's start with a triangular matrix. Remember
// that elements are stored in column-major order
// by default.
let t =
Matrix.CreateUpperTriangular(
4, 4,
[|
1.0; 0.0; 0.0; 0.0;
1.0; 2.0; 0.0; 0.0;
1.0; 4.0; 1.0; 0.0;
1.0; 3.0; 1.0; 2.0
|], MatrixElementOrder.ColumnMajor)
let b1 = Vector.Create(1.0, 3.0, 6.0, 3.0)
let b2 =
Matrix.CreateFromArray(4, 2,
[|
1.0; 3.0; 6.0; 3.0;
2.0; 3.0; 5.0; 8.0
|], MatrixElementOrder.ColumnMajor)
printfn "t = %O" (t.ToString("F4"))
//
// The Solve method
//
// The following solves m x = b1.
let x1 = t.Solve(b1)
printfn "x1 = %O" (x1.ToString("F4"))
// A second parameter specifies whether to overwrite
// the right-hand side with the result. Solve also
// returns the result, which we ignore here:
t.Solve(b1, true) |> ignore
printfn "b1 = %O" (b1.ToString("F4"))
// You can solve for multiple right hand side
// vectors by passing them in a DenseMatrix:
let x2 = t.Solve(b2, false)
printfn "x2 = %O" (x2.ToString("F4"))
//
// Related Methods
//
// You can verify whether a matrix is singular
// using the IsSingular method:
printfn "IsSingular(t) = %b" (t.IsSingular())
// The inverse matrix is returned by the Inverse
// method:
printfn "Inverse(t) = %O" (t.GetInverse().ToString("F4"))
// The determinant is also available:
printfn "Det(t) = %.4f" (t.GetDeterminant())
// The condition number is an estimate for the
// loss of precision in solving the equations
printfn "Cond(t) = %.4f" (t.EstimateConditionNumber())
printfn ""
//
// Symmetric systems and matrices
//
printfn "Symmetric matrices:"
// For the basics of working with symmetric
// matrices, see the SymmetricMatrices QuickStart
// Sample.
//
// Let's start with a symmetric matrix. Remember
// that elements are stored in column-major order
// by default.
let s =
Matrix.CreateSymmetric(4,
[|
1.0; 0.0; 0.0; 0.0;
1.0; 2.0; 0.0; 0.0;
1.0; 1.0; 2.0; 0.0;
1.0; 0.0; 1.0; 4.0
|], MatrixTriangle.Upper, MatrixElementOrder.ColumnMajor)
let b3 = Vector.Create(1.0, 3.0, 6.0, 3.0)
printfn "s = %O" (s.ToString("F4"))
//
// The Solve method
//
// The following solves m x = b1.
let x3 = s.Solve(b3)
printfn "x3 = %O" (x1.ToString("F4"))
// A second parameter specifies whether to overwrite the
// right-hand side with the result.
s.Solve(b2, true) |> ignore
printfn "b2 = %O" (b2.ToString("F4"))
// You can solve for multiple right hand side
// vectors by passing them in a DenseMatrix:
let x4 = s.Solve(b3, false)
printfn "x4 = %O" (x4.ToString("F4"))
//
// Related Methods
//
// You can verify whether a matrix is singular
// using the IsSingular method:
printfn "IsSingular(s) = %b" (s.IsSingular())
// The inverse matrix is returned by the Inverse
// method:
printfn "Inverse(s) = %O" (s.GetInverse().ToString("F4"))
// The determinant is also available:
printfn "Det(s) = %.4f" (s.GetDeterminant())
// The condition number is an estimate for the
// loss of precision in solving the equations
printfn "Cond(s) = %.4f" (s.EstimateConditionNumber())
printfn ""
//
// The CholeskyDecomposition class
//
// If the symmetric matrix is positive definite,
// you can use the CholeskyDecomposition class
// to optimize performance if multiple operations
// need to be performed. This class does the
// bulk of the calculations only once. This
// decomposes the matrix into G x transpose(G)
// where G is a lower triangular matrix.
//
// If the matrix is indefinite, you need to use
// the LUDecomposition class instead. See the
// LinearEquations QuickStart Sample for details.
printfn "Using Cholesky Decomposition:"
// The constructor takes an optional second argument
// indicating whether to overwrite the original
// matrix with its decomposition:
let cf = s.GetCholeskyDecomposition(false)
// The Factor method performs the actual
// factorization. It is called automatically
// if needed.
cf.Decompose()
// All methods mentioned earlier are still available:
let x5 = cf.Solve(b2, false)
printfn "x5 = %O" (x5.ToString("F4"))
printfn "IsSingular(m) = %b" (cf.IsSingular())
printfn "Inverse(m) = %O" (cf.GetInverse().ToString("F4"))
printfn "Det(m) = %.4f" (cf.GetDeterminant())
printfn "Cond(m) = %.4f" (cf.EstimateConditionNumber())
// In addition, you have access to the
// triangular matrix, G, of the composition.
printfn " G = %O" (cf.LowerTriangularFactor.ToString("F4"))
// Note that if the matrix is indefinite,
// the factorization will fail and throw a
// MatrixNotPositiveDefiniteException.
s.[0, 0] <- -99.0
let cf2 = s.GetCholeskyDecomposition()
try
cf2.Decompose()
with :? MatrixNotPositiveDefiniteException as e -> printfn "%s" e.Message
// To avoid this, you can use the TryDecompose method,
// which returns true if the decomposition succeeds,
// and false otherwise:
if (cf2.TryDecompose()) then
printfn "Decomposition succeeded!"
else
printfn "Decomposition failed!"
printf "Press Enter key to exit..."
Console.ReadLine() |> ignore