Skip to content

Commit

Permalink
Add LU and PLU decompositions. (#26)
Browse files Browse the repository at this point in the history 
* Add LU and PLU decompositions.

* Refactor & add description for LUP decomposition.

* Add docs for algorithms & fix some minor issues.

* Codestyle & format.

* Nullables added
  • Loading branch information
@TheSeems
TheSeems authored Jul 15, 2021
1 parent 8df6dda commit ebbdd1b
Show file tree
Hide file tree
Showing 13 changed files with 550 additions and 64 deletions.
3 changes: 2 additions & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,4 +3,5 @@
*.vs*
*/bin/*
*/obj/*
*/packages/*
*/packages/*
GenericTensor.sln.DotSettings.user
2 changes: 0 additions & 2 deletions GenericTensor.sln.DotSettings.user
View file

This file was deleted.

9 changes: 9 additions & 0 deletions GenericTensor/Core/Exceptions.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,4 +55,13 @@ public class InvalidDeterminantException : DataException
internal InvalidDeterminantException(string msg) : base(msg) {}
internal InvalidDeterminantException() : base() {}
}

/// <summary>
/// Thrown when there is no decomposition for provided matrix
/// </summary>
public class ImpossibleDecomposition : DataException
{
internal ImpossibleDecomposition(string msg) : base(msg) {}
internal ImpossibleDecomposition() : base() {}
}
}
139 changes: 115 additions & 24 deletions GenericTensor/Declaration.cs
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#region copyright

/*
* MIT License
*
Expand All @@ -22,6 +23,7 @@
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/

#endregion


Expand All @@ -40,6 +42,7 @@ public sealed partial class GenTensor<T, TWrapper>
: ICloneable where TWrapper : struct, IOperations<T>
{
#region Composition

/// <summary>
/// Creates a new axis that is put backward
/// and then sets all elements as children
Expand Down Expand Up @@ -67,6 +70,7 @@ public static GenTensor<T, TWrapper> Concat(GenTensor<T, TWrapper> a, GenTensor<
#endregion

#region Constructors

/// <summary>
/// Creates a tensor whose all matrices are identity matrices
/// <para>1 is achieved with <see cref="IOperations{T}.CreateOne"/></para>
Expand Down Expand Up @@ -145,36 +149,31 @@ public static GenTensor<T, TWrapper> CreateMatrix(int width, int height)
/// (its only argument is an array of integers which are indices of the tensor)
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> CreateTensor(TensorShape shape, Func<int[], T> operation, Threading threading = Threading.Single)
=> Constructors<T, TWrapper>.CreateTensor(shape, operation, threading);
public static GenTensor<T, TWrapper> CreateTensor(TensorShape shape, Func<int[], T> operation, Threading threading = Threading.Single) => Constructors<T, TWrapper>.CreateTensor(shape, operation, threading);

/// <summary>
/// Creates a tensor from an array
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> CreateTensor(T[] data)
=> Constructors<T, TWrapper>.CreateTensor(data);
public static GenTensor<T, TWrapper> CreateTensor(T[] data) => Constructors<T, TWrapper>.CreateTensor(data);

/// <summary>
/// Creates a tensor from a two-dimensional array
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> CreateTensor(T[,] data)
=> Constructors<T, TWrapper>.CreateTensor(data);
public static GenTensor<T, TWrapper> CreateTensor(T[,] data) => Constructors<T, TWrapper>.CreateTensor(data);

/// <summary>
/// Creates a tensor from a three-dimensional array
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> CreateTensor(T[,,] data)
=> Constructors<T, TWrapper>.CreateTensor(data);
public static GenTensor<T, TWrapper> CreateTensor(T[,,] data) => Constructors<T, TWrapper>.CreateTensor(data);

/// <summary>
/// Creates a tensor from an n-dimensional array
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> CreateTensor(Array data)
=> Constructors<T, TWrapper>.CreateTensor(data);
public static GenTensor<T, TWrapper> CreateTensor(Array data) => Constructors<T, TWrapper>.CreateTensor(data);

#endregion

Expand All @@ -186,6 +185,42 @@ public static GenTensor<T, TWrapper> CreateTensor(Array data)
/// </summary>
public GenTensor<T, TWrapper> RowEchelonFormSimple()
=> EchelonForm<T, TWrapper>.RowEchelonFormSimple(this);

/// <summary>
/// Decomposes a matrix into a triangular one.
/// Is of the Row Echelon Form (leading elements might be differ from ones).
///
/// In addition, returns a permutation that algorithm performs with rows.
/// Permutation array is size of rows there are in matrix.
///
/// Initial state of that array is:
/// 1 2 3 ... numberOfRows
///
/// For example, algorithm swaps first and third rows then:
/// 3 2 1 ... numberOfRows
///
/// It can be useful performing decompositions
/// </summary>
public (GenTensor<T, TWrapper>, int[]) RowEchelonFormPermuteSimple()
=> EchelonForm<T, TWrapper>.RowEchelonFormPermuteSimple(this);

/// <summary>
/// Decomposes a matrix into a triangular one.
/// Is of the Row Echelon Form (leading elements might be differ from ones).
///
/// In addition, returns a permutation that algorithm performs with rows.
/// Permutation array is size of rows there are in matrix.
///
/// Initial state of that array is:
/// 1 2 3 ... numberOfRows
///
/// For example, algorithm swaps first and third rows then:
/// 3 2 1 ... numberOfRows
///
/// It can be useful performing decompositions
/// </summary>
public (GenTensor<T, TWrapper>, int[]) RowEchelonFormPermuteSafeDivision()
=> EchelonForm<T, TWrapper>.RowEchelonFormPermuteSafeDivision(this);

/// <summary>
/// Decomposes a matrix into a triangular one.
Expand All @@ -196,7 +231,6 @@ public GenTensor<T, TWrapper> RowEchelonFormSafeDivision()
=> EchelonForm<T, TWrapper>.RowEchelonFormSafeDivision(this);



/// <summary>
/// Decomposes a matrix into a triangular one.
/// Is of the Row Echelon Form (leading elements are ones).
Expand All @@ -215,7 +249,6 @@ public GenTensor<T, TWrapper> RowEchelonFormLeadingOnesSafeDivision()
=> EchelonForm<T, TWrapper>.RowEchelonFormLeadingOnesSafeDivision(this);



/// <summary>
/// Finds the reduced echelon form of a matrix.
/// </summary>
Expand All @@ -228,6 +261,16 @@ public GenTensor<T, TWrapper> ReducedRowEchelonFormSimple()
/// </summary>
public GenTensor<T, TWrapper> ReducedRowEchelonFormSafeDivision()
=> EchelonForm<T, TWrapper>.ReducedRowEchelonFormSafeDivision(this);

/// <summary>
/// Finds the reduced echelon form of a matrix.
/// Uses safe division, i. e. perform division only when computing the final result.
///
/// Additionally returns row permutation
/// </summary>
public (GenTensor<T, TWrapper>, int[]) ReducedRowEchelonFormPermuteSafeDivision()
=> EchelonForm<T, TWrapper>.ReducedRowEchelonFormPermuteSafeDivision(this);

#endregion

#region Determinant
Expand Down Expand Up @@ -307,6 +350,7 @@ public void InvertMatrix()
/// </summary>
public void TensorMatrixInvert()
=> Inversion<T, TWrapper>.TensorMatrixInvert(this);

#endregion

#region Elementary matrix operations
Expand Down Expand Up @@ -398,70 +442,80 @@ public static GenTensor<T, TWrapper> TensorMatrixMultiply(GenTensor<T, TWrapper>
/// T1 + T2
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseAdd(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseAdd(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseAdd(a, b, threading);

/// <summary>
/// T1 - T2
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseSubtract(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseSubtract(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseSubtract(a, b, threading);

/// <summary>
/// T1 * T2
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseMultiply(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseMultiply(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseMultiply(a, b, threading);

/// <summary>
/// T1 / T2
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseDivide(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseDivide(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseDivide(a, b, threading);

/// <summary>
/// T1 + const
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseAdd(GenTensor<T, TWrapper> a, T b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseAdd(GenTensor<T, TWrapper> a, T b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseAdd(a, b, threading);

/// <summary>
/// T1 - const
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseSubtract(GenTensor<T, TWrapper> a, T b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseSubtract(GenTensor<T, TWrapper> a, T b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseSubtract(a, b, threading);

/// <summary>
/// const - T1
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseSubtract(T a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseSubtract(T a, GenTensor<T, TWrapper> b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseSubtract(a, b, threading);

/// <summary>
/// T1 * const
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseMultiply(GenTensor<T, TWrapper> a, T b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseMultiply(GenTensor<T, TWrapper> a, T b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseMultiply(a, b, threading);

/// <summary>
/// T1 / const
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseDivide(GenTensor<T, TWrapper> a, T b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseDivide(GenTensor<T, TWrapper> a, T b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseDivide(a, b, threading);

/// <summary>
/// const / T1
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GenTensor<T, TWrapper> PiecewiseDivide(T a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
public static GenTensor<T, TWrapper> PiecewiseDivide(T a, GenTensor<T, TWrapper> b,
Threading threading = Threading.Single)
=> PiecewiseArithmetics<T, TWrapper>.PiecewiseDivide(a, b, threading);

#endregion
Expand Down Expand Up @@ -493,6 +547,40 @@ public GenTensor<T, TWrapper> TensorMatrixPower(int power, Threading threading =

#endregion

#region Decompositions

/// <summary>
/// https://www.geeksforgeeks.org/l-u-decomposition-system-linear-equations/
/// </summary>
/// <returns>
/// LU decomposition
/// </returns>
public (GenTensor<T, TWrapper> l, GenTensor<T, TWrapper> u) LuDecomposition()
=> LuDecomposition<T, TWrapper>.Decompose(this);

/// <summary>
/// Find PLU decomposition: matrices P, L, U such that for original matrix A: PA = LU.
///
/// P stands for permutation matrix, permutations are made during the Gauss elimination process
/// L stands for LIBERTY lower triangle matrix
/// U stands for upper triangle matrix
///
/// Algorithm, given matrix A:
/// 1. Form an adjacent matrix (A|E)
/// 2. Find row echelon form of that matrix (U|L_0) and permutation of the rows
/// 3. Form permutation matrix P such that P_ij = \delta_{}
/// 4. Compute L = P * L_0^{-1}
///
/// Results are: P, L, U
/// </summary>
/// <returns>
/// LUP decomposition of given matrix
/// </returns>
public (GenTensor<T, TWrapper> p, GenTensor<T, TWrapper> l, GenTensor<T, TWrapper> u) PluDecomposition()
=> PluDecomposition<T, TWrapper>.Decompose(this);

#endregion

#region ToString & GetHashCode

/// <inheritdoc/>
Expand All @@ -517,7 +605,8 @@ public static GenTensor<T, TWrapper> VectorCrossProduct(GenTensor<T, TWrapper> a
/// <summary>
/// Calls VectorCrossProduct for every vector in the tensor
/// </summary>
public static GenTensor<T, TWrapper> TensorVectorCrossProduct(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b)
public static GenTensor<T, TWrapper> TensorVectorCrossProduct(GenTensor<T, TWrapper> a,
GenTensor<T, TWrapper> b)
=> VectorProduct<T, TWrapper>.TensorVectorCrossProduct(a, b);

/// <summary>
Expand Down Expand Up @@ -562,6 +651,7 @@ public GenTensor<T, TWrapper> Forward()
#endregion

#region Serialization

/*
* Serialization protocol:
*
Expand Down Expand Up @@ -595,6 +685,7 @@ public byte[] Serialize()
/// </returns>
public static GenTensor<T, TWrapper> Deserialize(byte[] data)
=> Serializer<T, TWrapper>.Deserialize(data);

#endregion
}
}
}
2 changes: 1 addition & 1 deletion GenericTensor/Functions/Determinant.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ internal static T DeterminantGaussianSafeDivision(GenTensor<T, TWrapper> t, int
return t.GetValueNoCheck(0, 0);

var n = diagLength;
var elemMatrix = EchelonForm<T, TWrapper>.InnerGaussianEliminationSafeDivision(t, n, n, out var swapCount);
var elemMatrix = EchelonForm<T, TWrapper>.InnerGaussianEliminationSafeDivision(t, n, n, null, out var swapCount);

var det = default(EchelonForm<T, TWrapper>.WrapperSafeDivisionWrapper<T, TWrapper>).CreateOne();
for (int i = 0; i < n; i++)
Expand Down
Loading

0 comments on commit ebbdd1b

Please sign in to comment.