help rowmat_utils

Title

aug

rowmat_utils -- Mata functions to manipulate (square) matrices in parallel collected as rows of a larger matrix

Contents

help entry Purpose

rm_matmult() Multiply (row) matrices

rm_matvecmult() Multiply (row) matrices and vectors

rm_vecvecmult() Multiply (row) vectors

rm_newtinv() Invert (row) matrices

rm_transpose() Transpose (row) matrices

rm_absrowsums() Compute (row) matrices' absolute value row sums

rm_abscolsums() Compute (row) matrices' absolute value column sums

rm_alpha0() Compute initial matrices for (newton) inversion

Description

The functions above are designed for use in situations in which one wishes to do a large number of
matrix operations in parallel. For example, suppose that one has a data set consisting of N observations,
where N is a large number, and, for each observation, matrix multiplication, transposition, inversion, etc. of some
n -by- n matrix is required. Rather than looping over N observations and inverting a matrix for each observation, the routines above might
instead be employed by rendering each nXn matrix as a 1-by-(n^2) row vector and applying the above operations.
The above-listed operations perform
the listed operations in parallel, looping over matrix entries, and therefore might run a bit faster than looping over observations.
For example, newtinv() could be used to compute inverses for a large
number of observations in parallel.

For examples as to how the above operations work, see the individual help entries.

Description of individidual operations

rm_matmult() -- Multiply in parallel square matrices with entries written as single rows

Syntax

real matrix rm_matmult(real matrix X, real matrix Y)

where

X: real matrix containing data (each row of which represents the entries, in order, of a square matrix)

Y: real matrix containing data (each row of which represents the entries, in order, of a square matrix)

Description

rm_matmult() multiplies sequences of matrices in parallel. The idea is, in situations in which a large number
of similar matrix multiplications must be carried out, to do so in parallel.

X and Y provide the sequences of matrices to be multiplied and must be of the same dimension.

Remarks

Examples:

Consider the sets of matrices:

. mata: X1=1,2,3 \ 4,5,6 \ 7,8,9
. mata: X2=2,3,4 \ 5,6,7 \ 8,9,10
. mata: Y1=3,4,5 \ 6,7,8 \ 9,10,11
. mata: Y2=4,5,6 \ 7,8,9 \ 10,11,12
. mata: X1*Y1
1 2 3
+
+
1 | 42 48 54 |
2 | 96 111 126 |
3 | 150 174 198 |
+
+

. mata: X2*Y2
1 2 3
+
+
1 | 69 78 87 |
2 | 132 150 168 |
3 | 195 222 249 |
+
+

. mata: X1=rowshape(X1,1)
. mata: X2=rowshape(X2,1)
. mata: X=X1\X2
. mata: Y1=rowshape(Y1,1)
. mata: Y2=rowshape(Y2,1)
. mata: Y=Y1\Y2
. mata: rm_matmult(X,Y)
1 2 3 4 5 6 7 8 9
+
+
1 | 42 48 54 96 111 126 150 174 198 |
2 | 69 78 87 132 150 168 195 222 249 |
+
+

Conformability

The rows of X and Y must be conformable square matrices. Thus, in general X,Y should
have 1,4,9,16,etc. rows and N columns.

Additional Comments

No special attempt at efficiently multiplying the matrices has been made;
like other rowmat_utils, rm_matmult()
might be useful in situations where the primary need is to do a large number of matrix multiplications and doing
them in parallel conveys an advantage over looping over matrices. That is, looping is over entries of the matrix, not over matrices.

rm_matvecmult() -- Multiply in parallel (square) matrices and vectors with entries written as single rows

Syntax

real matrix rm_matvecmult(real matrix X, real matrix y)

where

X: real matrix containing data (each row of which represents the entries, in order, of a square matrix)

y: real matrix containing data (each row of which represents the entries, in order, of a vector)

Description

rm_matvecmult() multiplies sequences of matrices and vectors in parallel. The idea is, in situations in which a large number
of similar matrix-vector multiplications must be carried out, to do so in parallel.

X and y provide the sequences of matrices and vectors, respectively,
to be multiplied and must be conformable.

Remarks

Examples:

Consider the sets of matrices and (column) vectors:

. mata: X1=1,2,3 \ 4,5,6 \ 7,8,9
. mata: X2=2,3,4 \ 5,6,7 \ 8,9,10
. mata: y1=3 \4 \ 5
. mata: y2=4\ 5 \ 6
. mata: X1*Y1
1
+------+
1 | 26 |
2 | 62 |
3 | 98 |
+------+

. mata: X2*Y2
1
+------+
1 | 47 |
2 | 92 |
3 | 137 |
+------+

. mata: X1=rowshape(X1,1)
. mata: X2=rowshape(X2,1)
. mata: X=X1\X2
. mata: y1=y1'
. mata: y2=y2'
. mata: Y=y1\y2
. mata: rm_matvecmult(X,y)
1 2 3
+
+
1 | 26 62 98 |
2 | 47 92 137 |
+
+

Conformability

X and y must be conformable; if the individual matrices represented by X is n-by-n (and therefore 1 -by- n^2
when written as a row, y should represent vectors that are n-by-1 (or 1 -by- n when
represented as a row). Thus, X should be of dimension N-by-n^2 and y should be
of dimension N-by-n.

Additional Comments

No special attempt at efficiently executing matrix-vector multiplication has
been attempted. Like other rowmat_utils, rm_matvecmult()
will be used in situations where the primary need is to do a large number of matrix-vector multiplications and doing
them in parallel conveys an advantage over looping. Looping is over entries of the matrix, not over matrices.

rm_vecvecmult() -- Scalar (dot) products of vectors in parallel, with entries written as single rows

Syntax

real matrix rm_vecvecmult(real matrix x, real matrix y)

where

x: real matrix containing data (each row of which represents the entries of a vector)

y: real matrix containing data (each row of which represents the entries of a vector)

Description

rm_vecvecmult() multiplies sequences of vectors in parallel. The idea is, in situations in which a large number
of similar vector-vector multiplications must be carried out, to do so in parallel.

x and y provide the sequences of vectors to be multiplied and must be conformable.

Remarks

Examples:

Consider the sets of matrices and (column) vectors:

. mata: x1=1 \2 \ 3
. mata: x2=2 \3 \ 4
. mata: y1=3 \4 \ 5
. mata: y2=4 \5 \ 6
. mata: x1'y1
26

. mata: x2'y2
47

. mata: x1=x1'
. mata: x2=x2'
. mata: x=x1\x2
. mata: y=y1\y2
. mata: rm_vecvecmult(x,y)
1
+----+
1 | 26 |
2 | 47 |
+----+

Conformability

x and y must be conformable (i. e., vectors of the same length). Thus, x should be of dimension N-by-n
and y should be of dimension N-by-n.

Additional Comments

No special attempt at efficiently executing vector-vector multiplication has
been attempted. Like other rowmat_utils, rm_vecvecmult()
will be useful in situations where the primary need is to do a large number of matrix-vector multiplications and doing
them in parallel conveys an advantage over looping. In fact, rm_vecvecmult() is identical to computing the row products
using the colon operator; in the example, this would be, rowsum(X:*Y).

rm_newtinv() -- Invert a sequence of matrices written as rows of a larger matrix in parallel using Newton iteration

Syntax

real matrix rm_newtinv(A, maxiter,
crit)

where

A: real matrix containing data (each row of which represents the entries, in order, of a square matrix)

maxiter: real scalar containing the number of maximum iterations rm_newtinv() will do in finding inverses. Often, a relatively
small number of iterations is sufficient (i.e., in the neighborhood of 20 maximum iterations)

crit: real scalar containing the convergence criterion; when improvements across all matrices are smaller than this value,
rm_newtinv() stops. A suitable criterion might be something in the neighborhoood of 1e-12.

Description

rm_newtinv() carries out simultaneous matrix inversion of a sequence of matrices which are all written as rows of a larger data matrix.
rm_newtinv() uses all the mata functions listed in rowmat_utils to accomplish this. The idea is to accomplish matrix inversion
of the sequence of matrices in parallel, avoiding having to loop over matrices to find inverses.

Remarks

The approach followed by rm_newtinv() is outlined in Pan and Schreiber (1990), who build on Shulz (1933), which originally proposed newton iteration
for computing the inverse of a nonsingular matrix A.

Consider a sequence of matrices, X[k],k=0,1,3,...,n. The inverse of A is sought. Newton inversion proceeds by iterating

X[k+1]=X[k](2I-AX[k])

Until convergence of the sequence of matrices X[k]. The resulting matrix X[k] is an (arbitrarily close) approximation of the inverse of A.
While there are more efficient and stable ways to invert a single matrix, Newton iteration has the advantage in that it is easily applied to many
matrices in parallel, and is relatively stable.

Of course, any iterative procedure requires a starting point, and convergence often depends upon good starting values. Pan and Schreiber (1990)
note that Ben-Israel and Cohen (1966) show that for a sufficiently small value of a scalar a0, X[0]=a0*A', results in convergence. While one can achieve more rapid
convergence with other choices of a, Pan and Schreiber (1990) show that choosing:

a0=1/(A_1*A_inf)

where A_1 denotes the largest of the absolute value of the row sums of the matrix, and A_inf is the largest of the absolute value of the column sums of
the matrix A (the p=1 norm and the p=infinity norm of the matrix A, respectively). The companion function rm_alpha0() calculates this starting value
for matrices written in row form.

Examples:

Consider the sets of matrices:

. mata: A1=2,5,1 \ 3,4,5 \ 4,4,2
. mata: A2=3,3,1 \ 1,2,3 \ 4,5,2
. mata: luinv(A1)
1 2 3
+
+
1 | -.2857142857 -.1428571429 .5 |
2 | .3333333333 0 -.1666666667 |
3 | -.0952380952 .2857142857 -.1666666667 |
+
+

. mata: luinv(A2)
1 2 3
+
+
1 | 1.833333333 .1666666667 -1.666666667 |
2 | -1.666666667 -.3333333333 1.333333333 |
3 | .5 .5 -.5 |
+
+

. mata: A1=rowshape(A1,1)
. mata: A2=rowshape(A2,1)
. mata: A=A1\A2
. mata: rm_newtinv(A,30,1e-12)
1 2 3 4 5 6 7 8 9
+
+
1 | -.2857142857 -.1428571429 .5 .3333333333 1.27202e-17 -.1666666667 -0.952380952 .2857142857 -.1666666667 |
2 | 1.833333333 .1666666667 -1.166666667 -1.666666667 -.333333333 1.333333333 .5 .5 .5 |
+
+

Conformability

The rows of A must be conformable square matrices. Thus, in general A should
have 1,4,9,16,etc. rows and N columns, where N is some number of matrices that need to be inverted.

Additional Comments

No special attempt at performing the necessary matrix operations has been made; if one is interested in how multiplication, transposition,
norm computation, etc. are carried out, see the entries in rowmat_utils. rm_newtinv() is not efficient for small jobs, but is intended for use
in situations where one must quickly invert a large number of relatively small matrices. Essentially rm_newtinv() replaces the task of looping over matrices
and inverting each one, with the combined tasks of successive iterations and looping over entries of the matrix.

rm_transpose() -- Transposes of an array of matrices written in single-row form

Syntax

real matrix rm_transpose(real matrix X)

where

X: real matrix containing data (each row of which represents the entries of a vector)

Description

rm_transpose() computes the transposes of a group of matrices, each of which is written in single-row form, in parallel. The
idea is to avoid looping over a large number of matrices.

Examples:

Consider the sets of matrices:

. mata: X1=1,2,3 \ 4,5,6 \ 7,8,9
. mata: X2=2,3,4 \ 5,6,7 \ 8,9,10
. mata: X1'
1 2 3
+
+
1 | 1 4 7 |
2 | 2 5 8 |
3 | 3 6 9 |
+
+
. mata: X2'
1 2 3
+
+
1 | 2 5 8 |
2 | 3 6 9 |
3 | 4 7 10 |
+
+
. mata: X1=rowshape(X1,1)
. mata: X2=rowshape(X2,1)
. mata: X=X1 \ X2
. mata: rm_transpose(X)
1 2 3 4 5 6 7 8 9
+
+
1 | 1 4 7 2 5 8 3 6 9 |
2 | 2 5 8 3 6 9 4 7 10 |
+
+

Conformability

X must be a list of square n-by-n matrices, which then occupy together a single N-by-n^2 matrix.

Additional Comments

The intent is that,
like other rowmat_utils, rm_transpose()
will be useful in situations where the primary need is to do a large number of matrix computations and doing
them in parallel conveys an advantage over looping. Essentially, rm_transpose() replaces the task of looping over matrices
with the task of looping over the entries of many matrices simultaneously.

rm_absrowsums() -- Row sums of absolute values of a sequence of matrices written in single-row form

Syntax

real matrix absrowsums(real matrix X)

where

X: real matrix containing data (each row of which represents the entries of a square matrix)

Description

rm_absrowsums() computes the row sums of the absolute values of a matrix, where each matrix is written in single-row form,
in parallel. The idea is to perform the operation while avoiding looping over a large number of matrices.

Examples:

Consider the sets of matrices and (column) vectors:

. mata: X1=1,-2,3 \ 4,5,-6 \ 7,-8,9
. mata: X2=2,-3,4 \ 5,-6,7 \ 8,-9,10
. mata: rowsum(abs(X1))
1
+----+
1 | 6 |
2 | 15 |
3 | 24 |
+----+
. mata: rowsum(abs(X2))
1
+----+
1 | 9 |
2 | 18 |
3 | 27 |
+----+

. mata: X1=rowshape(X1,1)
. mata: X2=rowshape(X2,1)
. mata: X=X1 \ X2
. mata: rm_absrowsums(X)
1 2 3
+
+
1 | 6 15 24 |
2 | 9 18 27 |
+
+

Conformability

X must be a list of square n-by-n matrices, which then occupy together a single N-by-n^2 matrix.

Additional Comments

Like other rowmat_utils, rm_absrowsums()
will be useful in situations where the primary need is to do a large number of matrix computations and doing
them in parallel conveys an advantage over looping. Looping is over entries in the matrix, not over matrices.

rm_abscolsums() -- Column sums of absolute values of a sequence of matrices written in single-row form

Syntax

real matrix rm_abscolsums(real matrix X)

where

X: real matrix containing data (each row of which represents the entries of a square matrix)

Description

rm_abscolsums() computes the row sums of the absolute values of a matrix, where each matrix is written in single-row form,
in parallel. The idea is to perform the operation while avoiding looping over a large number of matrices.

Examples:

Consider the sets of matrices and (column) vectors:

. mata: X1=1,-2,3 \ 4,5,-6 \ 7,-8,9
. mata: X2=2,-3,4 \ 5,-6,7 \ 8,-9,10
. mata: colsum(abs(X1))
1 2 3
+
+
1 | 12 15 18 |
+
+
. mata: colsum(abs(X2))
1 2 3
+
+
1 | 15 18 21 |
+
+

. mata: X1=rowshape(X1,1)
. mata: X2=rowshape(X2,1)
. mata: X=X1 \ X2
. mata: rm_abscolsums(X)
1 2 3
+
+
1 | 12 15 18 |
2 | 15 18 21 |
+
+

Conformability

X must be a list of square n-by-n matrices, which then occupy together a single N-by-n^2 matrix.

Additional Comments

Like other rowmat_utils, rm_abscolsums()
might be useful in situations where the primary need is to do a large number of matrix computations and doing
them in parallel conveys an advantage over looping. The operation is performed by looping over entries in the matrix, and not
looping over matrices.

rm_alpha0() -- Provides starting values for rm_newtinv().

For details and description see rm_newtinv() and rowmat_utils.

References

Ben-Israel, A. and D. Cohen. 1966. "On iterative computation of generalized inverses and associated projections," SIAM Journal on Numerical Analysis 3: 410-9.

Pan, V. and R. Schreiber. 1991. "An improved Newton iteration for the generalized inverse of a matrix, with applications," SIAM Journal of Scientific and Statistical Computing 12:1109-31.

Schultz, G. 1933. "Iterative berechnung der Reziproken Matrix," Zeitschrift fuer Angewandte Mathematik and Mechanik 13: 57-9.

Author

Matthew Baker, Hunter College and the Graduate Center, CUNY, mjbaker@hunter.cuny.edu