This section is devoted to very important practical techniques how to extract submatrices and how to build matrices from smaller matrices. Although all presented examples and codes are applied to matrices from space 𝔽^m×n with entries from field 𝔽 (which is either ℤ or ℚ or ℝ or ℂ), all these matrix operations can be extended for abstract matrices where entries could be arbitrary algebraic objects including lists.

https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/matrix
https://www.r-tutor.com/r-introduction/matrix

,p> Vectors play a fundamental roe in linear algebra since they are used to represent quantities of different algebraic structure. In R, vectors are essential data structures that facilitate various mathematical operations. This section will cover basic to advanced manipulations with vectors, providing clear explanations and R code examples.

Manipulations with vectors

Usually, an n-tuple (x₁, x₂, … , x_n), written in parenthesis and a row vector [x₁, x₂, … , x_n], written in brackets, look as the same object to human eyes. One of the pedagogical virtues of any software package is its requirement to pay close attention to the types of objects used in specific contexts. In particular, the computer algebra system Mathematica treats these two versions of a vector differently because

\[ \left( x_1 , x_2 , \ldots , x_n \right) \in \mathbb{F}^n , \]

but

\[ \left[ x_1 , x_2 , \ldots , x_n \right] \in \mathbb{F}^{n\times 1} , \]

where x [bold] is a matrix and 𝔽^m×n is read “m by n” indicating dimensions of a (rows by columns) matrix, not m multiplied by n. This denotes the space of m × n matrices with coefficients from field 𝔽.

To create a vector in R, use the c() function to combine elements.

You can find the length (number of elements) of a vector using the length() function

In R, you can create column vectors and row vectors using the matrix() function.

Similar operations lead to a row vector:

Example 1: You can convert a row vector into a column vector or vice versa using the t() function. This operation is called transposition.

■

End of Example 1

You can scale a vector by multiplying it by a scalar

Adding two vectors element-wise is straightforward in R.

A zero vector has all elements as zero, and a one vector has all elements as one.

Two vectors are orthogonal if their dot product is zero.

Normalization is the process of scaling a vector so that it has a length of 1. This is done by dividing each element of the vector by its magnitude (or norm).

You can extract specific elements or subvectors from a vector using indexing.

Now we present some examples with randomly generated entries.

Example 2: R provides functions to create vectors with randomly generated entities. Here’s an example using runif() to generate random numbers between 0 and 1.

Random complex numbers are also available Real and Imaginary may be separately accessed ■

End of Example 2

Another practical example:

Example 3: For three given vectors,

Clear[x, y, z]
x = {1, 2, 3};
y = {{1, 2, 3}};
z = {{1}, {2}, {3}};

Extracting entries from vectors:

x[[1]] (* extracting the first component *)

y[[1]]

{1, 2}

y[[1, 2]] (* extracting the first component *)

y[[1, 1 ;; 2]]

{1, 2}

z[[1]]

{1}

z[[2, 1]]

Let us randomly generate a 9-tuple:

x = RandomInteger[{0, 9}, 9]

{6, 8, 0, 5, 6, 2, 8, 3, 5}

Extracting a subvector with even indices may be created in several ways

x[[#]] & /@ Range[2, 9, 2]

{8, 5, 2, 3}

Select[PositionIndex[x][[All, 1]], EvenQ[#] &] // Keys

{8, 5, 2, 3}

Extract subvector with odd indices

x[[#]] & /@ Range[1, 9, 2]

6, 0, 6, 8, 5}

building subvector with even or odd entries and vector with third values as [1, 3, 6, 9, ..]

extract subvesto of the first left part

extract subvector of right part ■

End of Example 3

reverse order of entries in a matrix

Example 4:

■

End of Example 4

Example 5: Vectors can represent stock prices in a portfolio. You can calculate the overall value of a stock portfolio and analyze the contribution of each stock.

■

End of Example 5

Example 6: Vectors represent forces acting on an object. Calculate the net force acting on an object by summing up individual force vectors.

■

End of Example 6

We can visualize two dimensional vectors with R:

We can multiply vectors, but the output depends on the order of the vectors and their dimensions. There are known three (actually, there are more, but we consider only three the most popular ones) types of products. We start with the outer product of two vectors of arbitrary dimensions. Since the outer product of two vectors is a rectangular matrix, it does not matter to what space these vectors belong. However, it is convenient to consider one vector as a row, but another as a column.

The outer product (also known as tensor product) of two vectors x ∈ ℝ^m×1 and y ∈ ℝ^1×n is their matrix product whose entries are all products of an element in the first vector with an element from the second vector:

\[ {\bf x} \otimes {\bf y} = \left[ \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_m \end{array} \right] \,\begin{bmatrix} y_1 & y_2 & \cdots & y_n \end{bmatrix} = \begin{bmatrix} x_1 y_1 & x_1 y_2 & \cdots & x_1 y_n \\ x_2 y_1 & x_2 y_2 & \cdots & x_2 y_n \\ \vdots & \vdots & \ddots & \vdots \\ x_m y_1 & x_m y_2 & \cdots & x_m y_n \end{bmatrix} \in \mathbb{F}^{m\times n} . \]

The dot product of two vectors of the same size is the number

\[ \begin{bmatrix} x_1 & x_2 & \cdots & x_n \end{bmatrix} \bullet \begin{bmatrix} y_1 & y_2 & \cdots & y_n \end{bmatrix} = x_1 y_1 + \cdots + x_n y_n . \]

The dot product of two vectors is the sum of the products of their corresponding elements.

The Hadamard product of two vectors of the same size is the vector

\begin{align*} {\bf v} \odot {\bf u} &= \begin{bmatrix} v_1 & v_2 & \cdots & v_n \end{bmatrix} \odot \begin{bmatrix} u_1 & u_2 & \cdots & u_n \end{bmatrix} \\ &= \begin{bmatrix} v_1 u_1 & v_2 u_2 & \cdots & v_n u_n \end{bmatrix} . \end{align*}

Manipulations with Matrices

A matrices is usually defined in Mathematica using curly brackets:

m = {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}}

\( \displaystyle \quad \begin{pmatrix} 1&2&3&4 \\ 5&6&7&8 \\ 9&10&11&12 \end{pmatrix} \)

Extracting Blocks from Matrices

Example 7: Let us consider the 3-by-4 matrix

m = {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}}

\( \displaystyle \quad \begin{pmatrix} 1&2&3&4 \\ 5&6&7&8 \\ 9&10&11&12 \end{pmatrix} \)

At times we need just one row or column from a matrix

mRow2 = m[[2]]

{5, 6, 7, 8}

Extract row 3:

m[[3]]

Out[5]= {9, 10, 11, 12}

Now we perform operations with extracting columns:

mCol2 = m[[All, 2]]

{2, 6, 10}

Show column 2 in matrix form

m[[All, 2]] // MatrixForm

Out[6]= \( \displaystyle \quad \begin{pmatrix} 2 \\ 6 \\ 10 \end{pmatrix} \)

We can ask for a particular entry using the convention {r,c} to identify the row and column where the entry resides. Here is the second entry in the second row of our matrix

m[[2, 2]]

Going the other way, we can ask Mathematica to tell us where that element is in our matrix

Position[m, 6]

(2, 2)

We often need the elements in the diagonal of a matrix. First, we build a square matrix by using only the first three columns (matching the number of rows) of our sample matrix

m[[All, 1 ;; 3]]
Diagonal[%]

\( \displaystyle \quad \begin{pmatrix} 1&2&3 \\ 5&6&7 \\ 9&10&11 \end{pmatrix} \)

{1, 6, 11}

Notice this is the same as the diagonal of our original matrix because Mathematica starts at the upper left and ignores columns beyond what it needs for a Diagonal of the leftmost square portion of the matrix

Diagonal[m]

{1, 6, 11}

framedDiagonalMatrix[m_List] := Table[If[i == j, Framed[m[[i, j]], FrameStyle -> Red], m[[i, j]]], {i, Length[m]}, {j, Length[m[[1]]]}];
framedDiagonalMatrix[m]

\( \displaystyle \quad \begin{pmatrix} \fbox{1}&2&3&4 \\ 5&\fbox{6}&7&8 \\ 9&10&\fbox{11}&12 \end{pmatrix} \)

Replace a block of a matrix:

mat = RandomInteger[10, {5, 5}];
mat // MatrixForm

Out[2]= \( \displaystyle \quad \begin{pmatrix} 3&8&10&7&1 \\ 9&5&10&7&1 \\ 3&6&0&4&4 \\ 10&3&2&3&6 \\ 7&8&9&3&8 \end{pmatrix} \)

Update the 3-by-4 submatrix by using the short form of Span (;;) to specify the relevant span of rows and columns:

mat[[1 ;; 3, 1 ;; 4]] = {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};
mat // MatrixForm

Out[4]= \( \displaystyle \quad \begin{pmatrix} 0&0&0&0&1 \\ 0&0&0&0&1 \\ 0&0&0&0&4 \\ 10&3&2&3&6 \\ 7&8&9&3&8 \end{pmatrix} \)

Another version of Wolfram codes:

SeedRandom[39515];
mat = RandomInteger[10, {5, 5}]

\( \displaystyle \quad \begin{pmatrix} 2& 8& 7& 1& 4 \\ 0& 8& 2& 1& 3 \\ 9& 9& 0& 3& 7 \\ 2& 4& 9& 5& 0 \\ 0& 4& 8& 5& 9 \end{pmatrix} \)

Update the 3-by-4 submatrix by using the short form of Span (;;) to specify the relevant span of rows and columns:

mat2 = mat /. mat[[1 ;; 3, 1 ;; 4]] -> {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}

\( \displaystyle \quad \begin{pmatrix} 2&10&1&3&2 \\ 7&4&6&6&3 \\ 6&6&6&4&4 \\ 3&2&7&7&9 \\ 6&2&5&2&1 \end{pmatrix} \)

To pick out a submatrix, you can use Span (;;). First, define a 4×5 matrix:

mat = RandomInteger[10, {4, 5}];
mat // MatrixForm

Out[6]= \( \displaystyle \quad \begin{pmatrix} 0& 1& 7& 4& 5 \\ 1& 10& 4& 5& 4 \\ 4& 10& 1& 1& 8 \\ 5& 2& 4& 10& 2 \end{pmatrix} \)

SeedRandom[39326];
mat3 = RandomInteger[10, {5, 5}]

\( \displaystyle \quad \begin{pmatrix} 3& 1& 8& 2& 3 \\ 4& 0& 7& 5& 0 \\ 5& 5& 7& 7& 2 \\ 3& 10& 3& 0& 0 \\ 10& 6& 2& 10& 6 \end{pmatrix} \)

Extract the top-left 3×4 submatrix by using Span (;;) to specify the relevant span of rows and columns:

mat3[[1 ;; 3, 1 ;; 4]] // MatrixForm

\( \displaystyle \quad \begin{pmatrix} 3& 1& 8& 2 \\ 4& 0& 7& 5 \\ 5& 5& 7& 7 \end{pmatrix} \)

You can get the same block from a matrix with another command:

Take[mat3, {1, 3}, {1, 4}]

\( \displaystyle \quad \begin{pmatrix} 3& 1& 8& 2 \\ 4& 0& 7& 5 \\ 5& 5& 7& 7 \end{pmatrix} \)

Another option to omit or truncate some rows or/and columns. If one wants to delete, say the second row, type:

Drop[mat3, {2}]

\( \displaystyle \quad \begin{pmatrix} 3& 1& 8& 2& 3 \\ 5& 5& 7& 7& 2 \\ 3& 10& 3& 0& 0 \\ 10& 6& 2& 10& 6 \end{pmatrix} \)

Delete[mat3, 2]

\( \displaystyle \quad \begin{pmatrix} 3&1&8&2&3 \\ 5&5&7&7&2 \\ 3&10&3&0&0 \\ 10&6&2&10&6 \end{pmatrix} \)

If one wants to delete the first two rows, type:

Drop[mat3, 2]

\( \displaystyle \quad \begin{pmatrix} 5&5&7&7&2 \\ 3&10&3&0&0 \\ 10&6&2&10&6 \end{pmatrix} \)

Now suppose you want to remove certain rows and columns from a matrix. It make sense to write a short function:

remove[a_?MatrixQ, row_?VectorQ, col_?VectorQ] := Module[{nr, nc, krow, kcol}, {nr, nc} = Dimensions[a]; krow = Complement[Range[1, nr], row]; kcol = Complement[Range[1, nc], col]; a[[krow, kcol]] ]

As an example. consider a 4×5 matrix:

SeedRandom[39326];
b = RandomInteger[{0, 10}, {4, 5}];
MatrixForm[b]

Out[4]= \( \displaystyle \quad \begin{pmatrix} 3& 1& 8& 2& 3 \\ 4& 0& 7& 5& 0 \\ 5& 5& 7& 7& 2 \\ 3& 10& 3& 0& 0 \end{pmatrix} \)

Another Wolfram code:

SeedRandom[34676];
mat4 = RandomInteger[{0, 10}, {4, 5}]

\( \displaystyle \quad \begin{pmatrix} 2&7&10&10&1 \\ 0&2&10&7&2 \\ 1&9&7&5&2 \\ 10&3&7&4&9 \end{pmatrix} \)

Here we attempt to remove rows 2,5, and column 4. Note there is no row 5, so the command disregards it.

remove[b, {2, 5}, {4}] // MatrixForm

\( \displaystyle \quad \begin{pmatrix} 3& 1& 8& 3 \\ 5& 5& 7& 2 \\ 3& 10& 3& 0 \end{pmatrix} \)

Extract all elements except the outermost rows and columns (negative indices count from the end):

mat4[[2 ;; -2, 2 ;; -1]]

\( \displaystyle \quad \begin{pmatrix} 2&10&7&2 \\ 9&7&5&2 \end{pmatrix} \)

to be removed or modified:

============================================ finished

mat // MatrixForm
mat[[2 ;; -2, 2 ;; -1]]

Out[9]= \( \displaystyle \quad \begin{pmatrix} 3&1&8&5&0 \\ 0&0&8&6&8 \\ 8&0&0&10&0 \\ 4&4&8&9&8 \end{pmatrix} \)

Out[10]= {{0, 8, 6, 8}, {0, 0, 10, 0}}

Extract diagonal elements:

Diagonal[mat]

Out[11]= {3, 0, 0, 9}

Diagonal[mat, 2]

Out[12]= {8, 6, 0}

mat4

\( \displaystyle \quad \begin{pmatrix} 2&7&10&10&1 \\ 0&2&10&7&2 \\ 1&9&7&5&2 \\ 10&3&7&4&9 \end{pmatrix} \)

Diagonal[mat4]

{2, 2, 7, 4}

Note there are several diagonals, listed here from upper right to lower left

MatrixForm[#] &[Diagonal[mat4, #] & /@ {4, 3, 2, 1, 0, -1, -2, -3}]

\( \displaystyle \quad \left( \begin{array}{c} \{ 1 \} \\ \{ 10, \,2 \} \\ \{ 10, \,7, \, 2 \} \\ \{ 7, \, 10\, 5,\, 9 \} \\ \{ 2, \,2,\, 7,\, 4 \} \\ \{ 0,\, 9,\, 7 \} \\ \{ 1, \, 3 \} \\ \{ 10 \} \end{array} \right) \)

To find sum of diagonal elements (which is called trace), enter:

Total[Diagonal[mat]]

Out[13]= 12

Tr[mat4]

■

End of Example 7

Example 8: Let us consider the 3-by-4 matrix

m = {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}}

\( \displaystyle \quad \begin{pmatrix} 1&2&3&4 \\ 5&6&7&8 \\ 9&10&11&12 \end{pmatrix} \)

Adding two rows or columns. First, we add column 3 = column 3 + column 1:

m2 = m; m2[[All, 3]] += m2[[All, 1]];
m2 // MatrixForm

Out[8]= \( \displaystyle \quad \begin{pmatrix} 1&2&4&4 \\ 5&6&12&8 \\ 9&10&20&12 \end{pmatrix} \)

Now row 2 = row 2 + row 3:

m2 = m;
m2[[2]] += m2[[3]];
m2 // MatrixForm

Out[8]= \( \displaystyle \quad \begin{pmatrix} 1&2&3&4 \\ 14&16&18&20 \\ 9&10&11&12 \end{pmatrix} \)

The following picture is for fun and it is not related to the example.

MatrixPlot[{{1, 2, 3, 4}, {14, 16, 18, 20}, {9, 10, 11, 12}}]

Swapping rows or columns. Swap row 1 and row 3:

m2 = m; m2[[{1, 3}]] = m2[[{3, 1}]];
m2 // MatrixForm

Out[11]= \( \displaystyle \quad \begin{pmatrix} 9&10&11&12 \\ 5&6&7&8 \\ 1&2&3&4 \end{pmatrix} \)

Swap column 1 and 3:

m2[[All, {1, 3}]] = m2[[All, {3, 1}]];
m2 // MatrixForm

Out[13]= \( \displaystyle \quad \begin{pmatrix} 11&10&9&12 \\ 7&6&5&8 \\ 3&2&1&4 \end{pmatrix} \)

Multiply row 2 with 3:
To multiply row 2 of the original matrix by 3 multiply the matrix by a vector with 1 as the first and third element and 2 as the second element

m*{1, 3, 1} // MatrixForm

Out[14]= \( \displaystyle \quad \begin{pmatrix} 1&2&3&4 \\ 15&18&21&24 \\ 9&10&11&12 \end{pmatrix} \)

Multiply column 1 with 4:
Similarly, if you wish to multiply a single column by a scalar, leaving all other columns unchanged, here is column 1 of the original matrix times 4

((m // Transpose)*{4, 1, 1, 1}) // Transpose // MatrixForm

Out[15]= \( \displaystyle \quad \begin{pmatrix} 4&2&3&4 \\ 20&6&7&8 \\ 36&10&11&12 \end{pmatrix} \)

Swap row 1 and row 3:

m6 = m; m6[[{1, 3}]] = m6[[{3, 1}]]; Grid[{{"Old", "New"}, {MatrixForm[m], MatrixForm[m6]}}]

Old
\( \displaystyle \quad \begin{pmatrix} &&\mbox{old}& \\ 1&2&3&4 \\ 5&6&7&8 \\ 9&10&11&12 \end{pmatrix} \qquad \) \( \displaystyle \quad \begin{pmatrix} &&\mbox{New}& \\ 9&10&11&12 \\ 5&6&7 \\ 1&2&3&4 \end{pmatrix} \)

Swap column 1 and 3:

m7 = m6; m7[[All, {1, 2, 3}]] = m6[[All, {3, 2, 1}]]; Grid[{{"Old", "New"}, {MatrixForm[m6], MatrixForm[m7]}}]

Old
\( \displaystyle \quad \begin{pmatrix} &&\mbox{old}& \\ 9&10&11&12 \\ 5&6&7 \\ 1&2&3&4 \end{pmatrix} \qquad \) \( \displaystyle \quad \begin{pmatrix} &&\mbox{New}& \\ 11&10&9&12 \\ 7&6&5&8 \\ 3&2&1&4 \end{pmatrix} \)

■

End of Example 68

Building Matrices from Blocks

Example 7: Let us consider the 3-by-4 matrix

m = {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}}

\( \displaystyle \quad \begin{pmatrix} 1&2&3&4 \\ 5&6&7&8 \\ 9&10&11&12 \end{pmatrix} \)

Insert a column at position 2:

v = Range[21, 23]; Insert[m // Transpose, v, 2] // Transpose // MatrixForm

Out[4]= \( \displaystyle \quad \begin{pmatrix} 1& 21&2& 3& 4 \\ 5& 22&6& 7& 8 \\ 9& 23&10& 11& 12 \end{pmatrix} \)

Insert a row at position 2:

v = Range[30, 33]; Insert[m, v, 2] // MatrixForm

Out[8]= \( \displaystyle \quad \begin{pmatrix} 1&2&3&4 \\ 30&31&32&33 \\ 5&6&7&8 \\ 9&10&11&12 \end{pmatrix} \)

Alternately, we can insert a row at position 2

v2 = Range[30, 33]; m3 = Insert[m, v2, 2]

\( \displaystyle \quad \begin{pmatrix} 1&2&3&4 \\ 30&31&32&33 \\ 5&6&7&8 \\ 9&10&11&12 \end{pmatrix} \)

Position[m3, v2]

(2)

We insert a new column 3 = column 3 + column 1:

m4 = m;
m4[[All, 3]] += m4[[All, 1]];
Grid[{{"Old", "New"}, {MatrixForm[m], MatrixForm[m4]}}]

Now we insert into our original sample matrix a new row 2 = row 2 + row 3:

m5 = m;
m5[[2]] += m5[[3]];
Grid[{{"Old", "New"}, {MatrixForm[m], MatrixForm[m5]}}]

■

End of Example 7

Example 8: Let us consider the 3-by-4 matrix

m = {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}}

\( \displaystyle \quad \begin{pmatrix} 1&2&3&4 \\ 5&6&7&8 \\ 9&10&11&12 \end{pmatrix} \)

■

End of Example 8

A sparse matrix, or sparse array, is a matrix in which most elements are zero. A Band-Diagonal System generalizes the Tri-Diagonal condition, having non zero elements in a diagonal band and zero elements elsewhere.

SparseArray[{Band[{1, 1}] -> 1, Band[{1, 2}] -> 1, Band[{2, 1}] -> 1}, {5, 5}]

\( \displaystyle \quad \begin{pmatrix} 1&1&0&0&0 \\ 1&1&1&0&0 \\ 0&1&1&1&0 \\ 0&0&1&1&1 \\ 0&0&0&1&1 \end{pmatrix} \)

Example 9: Let us define a new sample matrix

SeedRandom[09654];
mat7 = Array[RandomInteger[{0, 100}] &, {4, 4}]

\( \displaystyle \quad \begin{pmatrix} 26&100&51&68& \\ 97&5&38&81 \\ 65&24&19&50 \\ 76&62&43&26 \end{pmatrix} \)

The matrix may be decomposed into its diagonal, upper trianglar and lower trianglar matrices

diag7 = DiagonalMatrix[Diagonal[mat7]];
ut7 = UpperTriangularMatrix[UpperTriangularize[mat7]];
lt7 = LowerTriangularMatrix[LowerTriangularize[mat7]];
Grid[{{"\nMatrix", "Diagonal\nMatrix", "Upper Triangular\nMatrix", "Lower Triangular\nMatrix"}, MatrixForm[#] & /@ {mat7, diag7, ut7, lt7}}, Alignment -> {Automatic, Center}]

Diagonal:
\( \displaystyle \quad \begin{pmatrix} 26&0&0&0 \\ 0&5&0&0 \\ 0&0&19&0 \\ 0&0&0&26 \end{pmatrix} \)

Upper Triangular:
\( \displaystyle \quad \begin{pmatrix} 26&100&51& 68 \\ 0&5&38&81 \\ 0&0&19&50 \\ 0&0&0&26 \end{pmatrix} \)

Lower Triangular:
\( \displaystyle \quad \begin{pmatrix} 26&0&0&0 \\ 97&5&0&0 \\ 65&24&19&0 \\ 76&62&43&26 \end{pmatrix} \)

Just as diagonals need not be in the absolute middle, triangular matrices need not be square and may be strictly or partially triangular, sometimes known as tri-diagonal

Grid[{{"Strictly Upper\nTrangular", "Upper\nTriangular", "\nTri-Diagonal"}, MatrixForm[#] & /@ {UpperTriangularize[mat7, 1], UpperTriangularize[mat7, 0], UpperTriangularize[mat7, -1]}}]

Strictly Upper Triangular:
\( \displaystyle \quad \begin{pmatrix} 0&100&51& 68 \\ 0&0&38&81 \\ 0&0&0&50 \\ 0&0&0&0 \end{pmatrix} \)

Upper Triangular:
\( \displaystyle \quad \begin{pmatrix} 26&100&51& 68 \\ 0&5&38&81 \\ 0&0&19&50 \\ 0&0&0&26 \end{pmatrix} \)

Tri-Diagonal:
\( \displaystyle \quad \begin{pmatrix} 26&100&51& 68 \\ 97&5&38&81 \\ 0&24&19&50 \\ 0&0&43&26 \end{pmatrix} \)

In linear algebra, the Crout matrix decomposition is an LU decomposition which decomposes a matrix into a lower triangular matrix (L), an upper triangular matrix (U), and, although not always needed, a permutation matrix (P). It was developed by the American mathematician Prescott Durand Crout (1907--1984). The Crout matrix decomposition algorithm differs slightly from the Doolittle method. Doolittle's method returns a unit lower triangular matrix and an upper triangular matrix, while the Crout method returns a lower triangular matrix and a unit upper triangular matrix.

crouts[matrix_] := Module[{n = Length[matrix], L, U}, L = IdentityMatrix[n];
U = DiagonalMatrix[Diagonal[matrix]];
Do[ L[[i, j]] = matrix[[i, j]] - Sum[L[[i, k]]*U[[k, j]], {k, 1, j - 1}]; U[[j, i]] = (matrix[[j, i]] - Sum[L[[j, k]]*U[[k, i]], {k, 1, j - 1}])/L[[j, j]], {j, 1, n}, {i, j + 1, n} ]; {L, U} ];

We demonstrate applications of Crout's algorithm:

cr7 = crouts[mat7];
Grid[{{"Lower Triangular\nMatrix", "Upper Triangular\nMatrix"}, MatrixForm[#] & /@ (cr7[[#]] & /@ Range[2])}, Alignment -> {Automatic, Center}]

Lower Triangular:
\( \displaystyle \quad \begin{pmatrix} 1&0&0&0 \\ 97&1&0&0 \\ 65&6476&1&0 \\ 76&7538&-37007875 & 1 \end{pmatrix} \)

Upper Triangular:
\( \displaystyle \quad \begin{pmatrix} 26&100&51&68 \\ 0&5&-4909 &6515 \\ 0&0&19&-42195510 \\ 0&0&0&26 \end{pmatrix} \)

Roger check
Mathematica helps with verification:

lower = {{1, 0, 0, 0}, {97, 1, 0, 0}, {65, -6476, 1, 0}, {76, -7538, -37007875, 1}};
upper = {{26, 100, 51, 68}, {0, 5, -4909, -6515}, {0, 0, 19, -42195510}, {0, 0, 0, 26}};
lower.upper

\( \displaystyle \quad \begin{pmatrix} 26&100&51&68 \\ 2522&9705&38&81 \\ 1690& -25880& 31794018 \\ 1976& -30090& -666141707& 1561566208756514 \end{pmatrix} \)

upper.lower

\( \displaystyle \quad \begin{pmatrix} 18209& -842760& -2516535449& 68 \\ -813740& 80900759& 241106300716& -6515 \\ -3206857525& 318069631336& 1561566159641269& -42195510 \\ 1976& -195988& -962204750& 26 \end{pmatrix} \)

The built-in Mathematic function LUDecomposition interchanges rows prior to decomposition to provide numerical stability; we demonstrate this function in a special section. ■

End of Example 9

Query

Axier, S., Linear Algebra Done Right. Undergraduate Texts in Mathematics (3rd ed.). Springer. 2015, ISBN 978-3-319-11079-0.
Beezer, R.A., A First Course in Linear Algebra, 2017.

Introduction to Linear Algebra

Systems of Linear Equations

Matrix Algebra

Vector Spaces

Eigenvalues, Eigenvectors

Euclidean Spaces

Matrix Decompositions

Applications

Functions of Matrices

Miscellany

Preliminaries

Glossary

Reference

Manipulations with vectors

Manipulations with Matrices

Extracting Blocks from Matrices

Building Matrices from Blocks

Query