An Orthogonally Based Pivoting Transformation of Matrices and Some Applications

In this paper we discuss the power of a pivoting transformation introduced by Castillo, Cobo, Jubete, and Pruneda [Orthogonal Sets and Polar Methods in Linear Algebra: Applications to Matrix Calculations, Systems of Equations and Inequalities, and Linear Programming, John Wiley, New York, 1999] and its multiple applications. The meaning of each sequential tableau appearing during the pivoting process is interpreted. It is shown that each tableau of the process corresponds to the inverse of a row modified matrix and contains the generators of the linear subspace orthogonal to a set of vectors and its complement. This transformation, which is based on the orthogonality concept, allows us to solve many problems of linear algebra, such as calculating the inverse and the determinant of a matrix, updating the inverse or the determinant of a matrix after changing a row (column), determining the rank of a matrix, determining whether or not a set of vectors is linearly independent, obtaining the intersection of two linear subspaces, solving systems of linear equations, etc. When the process is applied to inverting a matrix and calculating its determinant, not only is the inverse of the final matrix obtained, but also the inverses and the determinants of all its block main diagonal matrices, all without extra computations.


Introduction.
Castillo, Cobo, Fernández-Canteli, Jubete, and Pruneda [2] and Castillo, Cobo, Jubete, and Pruneda [3] have recently introduced a pivoting transformation of a matrix that has important properties and has been shown to be very useful to solve a long list of problems in linear algebra. The aim of this paper is to show the power of this transformation, clarify the meaning of the partial results obtained during the computationl process, and illustrate the wide range of applications of this transformation to solve common problems in linear algebra, such as calculating inverses of matrices, determinants or ranks, solving systems of linear equations, etc.
The reader interested in a classical treatment of these problems can, for example, consult the works of Burden and Faires [1], Golub and Van Loan [5], Gill et al. [6], and Press et al. [8].
The new methods arising from this transformation have complexity identical to that associated with the Gauss elimination method (see Castillo, Cobo, Jubete, and Pruneda [3]). However, they are specially suitable for updating solutions when changes in rows, columns, or variables are done. In fact, when changing a row, column, or variable, a single step of the process allows one to obtain (update) the new solution without the need to start again from scratch. For example, updating the inverse of an n × n matrix when a row is changed requires one instead of n steps, a drastic reduction in computational power.
In this paper we introduce the pivoting transformation and its applications only from the algebraic point of view. Discussing the numerical properties and performance of this method with respect to stability, ill conditioning, etc., which must be done carefully and taking into account its applications (see Demmel [4] and Higham [7]), will be the aim of another paper.
The paper is structured as follows. In section 2 the pivoting transformation is introduced. In section 3 its main properties are discussed. In section 4 an orthogonalization algorithm is derived. In section 5 some applications are given and illustrated with examples. Finally, some conclusions are given in section 6.

Pivoting transformation.
The main tool to be used in this paper consists of the so-called pivoting transformation, which transforms a set of vectors V j = {v j 1 , . . . , v j n } into another set of vectors V j+1 = {v j+1 where t j j = 0 and t j k , k = j, are arbitrary real numbers. In what follows we consider that the vectors above are the columns of a matrix V j .
This transformation can be formulated in matrix form as follows. Given a matrix where v j i , i = 1, . . . , n, are column vectors, a new matrix V j+1 is defined via where M −1 j is the inverse of the matrix M j = (e 1 , . . . , e j−1 , t j , e j+1 , . . . , e n ) T , (2.3) where e i is the ith column of the identity matrix, the transpose of t j being defined by for some predestined vector u j .
Since t j j = 0, the matrix M j is invertible. It can be proved that M −1 j is the identity matrix with its jth row replaced by This transformation is used in well-known methods, such as the Gaussian elimination method. However, different selections of the t-values lead to completely different results. In this paper we base this selection on the concept of orthogonality, and assume a sequence of m transformations associated with a set of vectors {u 1 , . . . , u m }. Downloaded 04/23/13 to 193.144.185.29. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 3. Main properties of the pivoting transformation. As we shall see, the pivoting transformation has very important and useful properties that are illustrated in the following theorems.
The first theorem proves that given a matrix V, the pivoting transformation transforms its columns without changing the linear subspace they generate.
implies the other way. In fact, this theorem is true by merely looking at (2.2).
The following theorem shows that the pivoting process (2.2) with the pivoting strategy (2.4) leads to the orthogonal decomposition of the linear subspace generated by the columns of V j with respect to vector u.
Theorem 3.2 (orthogonal decomposition with respect to a given vector). Assume now a vector u j = 0 and let t j In addition, the linear subspace orthogonal to and its complement is L(v j+1 j ). In other words, the transformation (2.2) gives the generators of the linear subspace orthogonal to u j and the generators of its complement.
Proof. This theorem follows quickly from (2.4) and (2.2) because Finally, Theorem 3.1 guarantees that L(v j+1 j ) is the complement. Remark 1. Note that Theorem 3.2 allows us to obtain the linear subspace orthogonal to a given vector u j in any case. If t j j = 0, we can reorder the v vectors until we satisfy the condition t j j = 0 or we find that t j j = 0 ∀j = 1, . . . , n, in which case the orthogonal set to The following two theorems show that the pivoting transformation (2.2) allows obtaining the linear space orthogonal to a given linear space, in another linear space. Theorem 3.3. If we sequentially apply the transformation in Theorem 3.1 based on a set of linearly independent vectors {u 1 , . . . , u j }, the orthogonalization and normalization properties in (3.1) are kept. In other words, we have where δ rk are the Kronecker deltas.
Proof. We prove this by induction over j.
Step 1. The theorem is true for j = 1, because from (3.1) we have u T 1 v 2 k = δ 1k . Downloaded 04/23/13 to 193.144.185.29. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Step j. We assume that the theorem is true for j, that is, Step j + 1. We prove that it is true for j + 1. In fact, we have Theorem 3.4 (orthogonal decomposition with respect to a given linear subspace). Assume the linear subspace L{u 1 , . . . , u n }. We can sequentially use Theorem 3.2 to obtain the orthogonal set to L{u 1 , . . . , u n } in a given subspace L(V 1 ). Let t j i be the dot product of u j and v j i . Then assuming, without loss of generality, that t q−1 q = 0, we obtain In addition, we have The proof can easily be obtained using Theorem 3.3.
The following remarks point out the practical significance of the above four theorems.
Remark 2. The linear subspace orthogonal to the linear subspace generated by vector u j is the linear space generated by the columns of V k for any k ≥ j + 1 with the exception of its pivot column, and its complement is the linear space generated by this pivot column of V k for any k ≥ j + 1.
Remark 3. The linear subspace, in the linear subspace generated by the columns of V 1 , orthogonal to the linear subspace generated by any subset W = {u k |k ∈ K} is the linear subspace generated by the columns of V , ≥ max k∈K k + 1, with the exception of all pivot columns associated with the vectors in W , and its complement is the linear subspace generated by the columns of V , ≥ max k∈K k + 1, which are their pivot columns.

The orthogonalization algorithm.
In this section we describe an algorithm for obtaining orthogonal decompositions, which is based on Theorem 3.4. Algorithm 1.
Step 3: Calculate the dot products t j = u T i w j , j = 1, . . . , s. Downloaded 04/23/13 to 193.144.185.29. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Step 4: For j = +1 to s locate the pivot column r as the first column not orthogonal to u i , that is, t r = 0. If there is no such a column go to Step 7. Otherwise, continue with Step 5.
Step 5: Increase in one unit, divide the r column by t r , and if r = , switch columns and r and associated dot products t r and t .
Step 6: For j = 1 to s and j = r do the following: Step 7: If i = m, go to Step 8. Otherwise, increase i in one unit and go to Step 3.
Step 8: Return L(W 2 ) = L (w , . . . , w s ) as the orthogonal subspace of L(U) in L(V 1 ) and L(W 1 ) = L (w 1 , . . . , w −1 ) as its complement. Remark 4. If the pivoting process were used taking into account numerical considerations, Step 4 should be adequately modified by the corresponding pivot selecting strategy (maximum pivot strategy, for example). In this case, the corresponding permutation in Step 8 is required. Note that in this paper only algebraic considerations are used.
The process described in Algorithm 1 can be organized in a tabular form. A detailed description is given in the following example.
Example 1 (orthogonal decomposition). Consider the linear subspace of L(V 1 ) = R 5 : We organize the procedure in a tabular form (see Table 4.1). First, to obtain the orthogonal decomposition of L(V 1 ) with respect to L(U), we construct the initial tableau (see Iteration 1 in Table 4.1), starting with the identity matrix V 1 . The first column of this table is the first generator of L(U) and the generators of the subspace to be decomposed are in the other columns. The last row contains the inner products of the vector in the first column by the corresponding column vectors.
Next, the first nonnull element in the last row is identified and the corresponding column is selected as the pivot column, which is boldfaced in Iteration 1. Downloaded 04/23/13 to 193.144.185.29. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Finally, it is necessary to perform the pivoting process and to update the first column and the last row of the table with the next generator of L(U) and the new inner products. Then, we get the second table (see Iteration 2 in Table 4.1). In order to select the pivot column, we have to look for the first nonnull element in the last row starting with its second element because we are in the second iteration. Then, the selected column is the third one, and before performing the pivoting process, interchange of second and third columns must be done.
We repeat the pivoting process, incorporate the last generators of L(U), and obtain the new dot products. We select the pivot column, starting at column three, and look for a nonnull dot product, obtaining the fourth column as the pivot. Next, we repeat the normalization and pivoting processes and, finally, we get the Output tableau in Table 4.1, where the first three vectors are the generators of the complement subspace and the last two are the generators of the orthogonal subspace. Italicized columns are used in all iterations to refer to the complementary subspace.
Thus, the orthogonal decomposition becomes Note that, from the Output tableau, we can obtain the linear subspace orthogonal to the linear subspace generated by any subset of the initial set of vectors. For example, the orthogonal complement of the linear subspace generated by the set Table 4.1) which can also be written as (see Iteration 3 in Table 4.1) Similarly,

Applications.
In addition to obtaining the linear subspace orthogonal to a linear space generated by one or several vectors in a given linear subspace, the proposed orthogonal pivoting transformation allows solving the following problems: 1. calculating the inverse of a matrix, 2. updating the inverse of a matrix after changing a row, 3. determining the rank of a matrix, 4. calculating the determinant of a matrix, 5. updating the determinant of a matrix after changing a row, 6. determining whether or not a set of vectors is linearly independent, 7. obtaining the intersection of two linear subspaces, 8. solving a homogeneous system of linear equations, 9. solving a complete system of linear equations, 10. deciding whether or not a linear system of equations is compatible.

Calculating the inverse of a matrix.
The following theorem shows that Algorithm 1 can be used for obtaining the inverse of a matrix.
Theorem 5.1. Assume that Algorithm 1 is applied to the rows of matrix A = (a 1 , . . . , a n ) T using a nonsingular initial matrix V 1 . Then the matrix whose columns Downloaded 04/23/13 to 193.144.185.29. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php are in the last tableau V n+1 is the inverse of matrix A. In addition, if we start with V 1 being the identity matrix, in the process we obtain the inverses of all block main diagonal matrices.
Proof. Matrices V j for j = 2, . . . , n + 1 are obtained, using the transformations This proves that A −1 = V n ; that is, the inverse of A is the matrix whose columns are in the final tableau obtained using Algorithm 1.
The second part of the theorem is obvious because the lower triangular part of the identity matrix is null and does not affect the dot products and the pivoting transformations involved in the process.
Example 2 (matrix inverses). Consider the following matrix A, where the block main diagonal matrices are shown, and its inverse A −1 : Table 5.1 shows the iterations for inverting A using Algorithm 1. The inverse matrix A −1 is obtained in the last iteration (see Table 5.1). In addition, Table 5.1 also contains the inverses of the block main diagonal matrices indicated below (see the marked matrices in Iterations 2 to 5 in Table 5.1). The important result is that this is obtained with no extra computation.
Finally, we mention that the 5 × 5 matrices we obtain in Iterations 2 to 5 in Table  5.1 are the inverses of the matrices that result from replacing in the unit matrix its rows by the rows of matrix A. For example, the matrix in Iteration 4 is such that Example 3 (inversion of a matrix starting from a regular matrix). The proposed pivoting process can be done starting with an arbitrary nonsingular matrix. For example, if we start with the matrix we get the results in Table 5.2, i.e., the same inverse. Downloaded 04/23/13 to 193.144.185.29. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Updating the inverse of a matrix after changing a row.
In this section we start by giving an interpretation to each tableau obtained in the inversion process of a matrix.
Since, according to Theorem 3.3, the pivoting transformation does not alter the orthogonal properties of previous vectors, we can update the inverse of a matrix after changing a row by an additional pivoting transformation in which the new row vector is used.
To illustrate this, we use the results in Table 5.2. Note that the matrices in Iterations 2 to 5 correspond to the matrices obtained from B −1 after sequentially replacing the row which number coincides with the number of the pivot column by their associated u-vectors. In other words, matrices in Table 5.2, Iterations 2 to 5, are the inverses of the following matrices:

Determining the rank of a matrix.
In this section we see that Algorithm 1 also allows one to determine the rank of a matrix. Downloaded 04/23/13 to 193.144.185.29. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php  In an n-dimensional linear space, the rank of a matrix U coincides with n minus the dimension of its orthogonal complement. Thus, if during the orthogonalization process we start with a nonsingular matrix as columns and we can find a pivot in all iterations, then the corresponding matrix is full rank. Otherwise, the rank is equal to the number of pivot columns we can find.
Example 4 (rank of a matrix). Assume that we are interested in calculating the rank of the matrix In Table 5.3 we show the iterations for obtaining its rank. We can see that the rank of A is 3, since the third and fifth iterations have no pivots.

Calculating the determinant of a matrix.
The following theorem shows that the determinant of a matrix can be calculated by means of Algorithm 1.
Theorem 5.2 (determinant of a matrix). The determinant of a matrix A can be calculated by Algorithm 1 by multiplying the normalizing constants t r , l = 1, . . . , n, used in Step 5 and (−1) p , where p is the number of interchanges of columns that have occurred when executing the algorithm. If we start the algorithm with the identity matrix, the determinants of the block main diagonal matrices referred to in section 5.1 are the partial products.
Proof. Assume that we start in Step 1 with an identity matrix, W = I n , which has a determinant of one. In the inverting process, we transform this matrix using two different transformations: the pivoting step, which does not alter the determinant Downloaded 04/23/13 to 193.144.185.29. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php value of the matrix, and the normalization step, which divides its determinant by t r (see (2.1)). In addition, we multiply it by −1 each time we switch columns. Since If we start with an identity matrix, the lower triangular part of the identity matrix is null and does not affect the dot products and the pivoting transformations involved in the process. Thus, the result holds.
Example 5 (determinant of a matrix). The determinant of the matrix in Example 2 is obtained by multiplying the normalizing constants, that is, the last values in the boldfaced columns in Table 5.1. Thus, we have The determinants of the block main diagonal matrices in (5.3) are 1, 2, 2, 2, and 7, respectively.
Remark 5. If instead of starting with the identity matrix I n , we start with a nonsingular matrix B with determinant |B|, expression (5.6) becomes

Updating the determinant of a matrix after changing a row.
According to (5.6) or (5.7), the determinant is updated by multiplying the previous determinant by the dot product of the new row by the associated pivot column.
Example 6 (updating the determinant after changing a row). Consider the matrix A and its inverse A −1 ,  and assume that we want to calculate the determinant of the matrix Since |A| = 7, we have

Determining whether or not a set of vectors is linearly dependent.
To know whether or not a set of vectors is linearly independent, we use the property Thus, the problem reduces to obtaining a set of generators of the orthogonal complement of L{u 1 , . . . , u n−1 } and checking that the dot products of each of its generators by u n are null.
Note that this problem is the same as determining whether or not a vector belongs to a linear subspace. If we use the pivoting process (see Table 5.4), we have no problem finding a pivot column for the first three vectors, but there is no pivot column for the fourth vector. This means that the fourth vector is a linear combination of the first three.

5.7.
Obtaining the intersection of two linear subspaces. Theorem 3.4 allows us to obtain the intersection of two linear subspaces S 1 and S 2 by noting that In fact, we can obtain first S ⊥ 2 , the orthogonal to S 2 , by letting L(V 1 ) = R n in Theorem 3.4 and then find the orthogonal to S ⊥ 2 in S 1 , using S 1 as L(V 1 ). Alternatively, we can obtain first S ⊥ 1 , the orthogonal to S 1 , by letting L(V 1 ) = R n in Theorem 3.4 and then find the orthogonal to S ⊥ 1 in S 2 , using S 2 as L(V 1 ). Downloaded 04/23/13 to 193.144.185.29. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Example 8 (intersection of moving subspaces). Consider the linear subspaces and assume that we wish to 1. determine the intersection Q 1 = S 1 ∩ S 2 for all values of the time parameter 0 ≤ t ≤ 2π; 2. find the t-values for which we have Then we have the following: 1. By definition we can write Using the procedure in Theorem 3.4 and starting with v 1 , we get (5.10) Since p = 0 ∀t, using the orthogonalization procedure in Theorem 3.4, we obtain and proceeding with v 4 and taking into account that

By a similar process we get
Since Q 2 is orthogonal to v 2 and Q 1 = Q 2 , then Q 1 is orthogonal to v 2 and, in particular, u 2 is orthogonal to v 2 . Similarly, since Q 1 is orthogonal to v 4 and Q 1 = Q 2 , then Q 2 is orthogonal to v 4 and, in particular, u 3 is orthogonal to v 4 ; that is, and, conversely, Thus, we get
Example 10 (a complete system of linear equations). Consider the system of equations which, using the auxiliary variable x 5 , can be written as (5.13). Since the solution of the homogeneous system (5.13) was already obtained, now we only need to force x 5 = 1 and return to the initial set of variables. Thus, the solution is (x 1 , x 2 , x 3 , x 4 ) = (0, 2, 0, 0) + ρ 1 (1, −1, 1, 1) , (5.20) where ρ 1 is an arbitrary real number.
Assume that now we add to system (5.19) the equation The new solution is obtained by an extra pivoting step using the new vector (0, 1, 0, −1, 0) (see Table  5.6).
Thus, analyzing the compatibility of the system of equations (5.15) reduces to finding the orthogonal complement L {w 1 , . . . , w p } of L {a 1 , . . . , a n } and checking whether or not bW T = 0.
The computational process arising from this procedure has a complexity which coincides exactly with that for the Gauss elimination procedure. However, it has one important advantage: W is independent of b and so we can analyze the compatibility of any symbolic vector u n+1 without extra computations.
Example 11 (compatibility of a linear system of equations). Suppose that we are interested in determining the conditions under which the system of equations is compatible. Then, using Algorithm 1, we get (see Table 5

Conclusions.
A pivoting transformation, based on the orthogonality concept, has been discussed and some of its applications to solve common linear algebra problems have been given. The main advantage of the suggested method with respect to the Gauss elimination method is that the intermediate results arising in the solution process are easily interpretable. This leads to immediate methods to update solutions of several problems, such as calculating the inverse or the determinant of a matrix, solving system of linear equations, etc., when small changes are done (changes in rows, columns, and/or variables). When the method is applied to inverting a matrix and calculating its determinant, not only is the inverse of the final matrix obtained, but also the inverses and determinants of all its block main diagonal matrices, without extra computations.