[Numpy-discussion] Implementing elementary matrices

[Matt Pagan]

Greetings!
I made a patch for NumPy that adds a function for easily creating
elementary matrices. Sorry for not knowing the process for submitting
patches.

Is this function something the NumPy community could see adding to the
codebase? Are there ways I can improve on this?

diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py
index 12c0f9b..10073af 100644
--- a/numpy/lib/twodim_base.py
+++ b/numpy/lib/twodim_base.py
@@ -967,3 +967,85 @@ def triu_indices_from(arr, k=0):
     if not (arr.ndim == 2 and arr.shape[0] == arr.shape[1]):
         raise ValueError("input array must be 2-d and square")
     return triu_indices(arr.shape[0], k)
+
+def elem(N, i, j=None, t=None, dtype=float):
+    """
+    Return an elementary matrix.
+
+    Parameters
+    ----------
+    N : int
+        The size of the NxN array to be returned. Elementary matrices
+        should be square.
+    i : int
+        The index of the first row on which operations are to be
+        performed.
+    j : int
+        If set, the index of the second row of which operations are to
+        be performed.
+    t : scalar
+        If set, the factor by which a given row will be multiplied.
+
+    Returns
+    -------
+    m: ndarray of shape (NxN)
+       The identity matrix after a single row operation has been
+       performed on it.
+
+    See also
+    --------
+    eye, identity
+
+    Examples
+    -------
+    To swap the the first and third rows of a 4x4 identity matirx:
+
+    >>> L = elem(4, 0, 2)
+    >>> L
+    array([[ 0.,  0.,  1.,  0.],
+           [ 0.,  1.,  0.,  0.],
+           [ 1.,  0.,  0.,  0.],
+           [ 0.,  0.,  0.,  1.]])
+
+    This array then becomes quite useful for matrix multiplication.
+
+    >>> H = np.matrix([[ 2,  3,  5,  7],
+           [11, 13, 17, 19],
+           [23, 29, 31, 37],
+           [41, 43, 47, 53]])
+    >>> L*H
+    matrix([[ 23.,  29.,  31.,  37.],
+            [ 11.,  13.,  17.,  19.],
+            [  2.,   3.,   5.,   7.],
+            [ 41.,  43.,  47.,  53.]])
+
+    When the elemntary matrix is multiplied by the given matrix, the
+    result is the given matrix with it's first and third rows swapped.
+
+    If the given matrix is multiplied by the elementary matrix (i.e.,
+    the multiplication takes place in reverse order, the result is
+   the given matrix with its first and third columns swapped.
+
+    >>> H*L
+    matrix([[  5.,   3.,   2.,   7.],
+            [ 17.,  13.,  11.,  19.],
+            [ 31.,  29.,  23.,  37.],
+            [ 47.,  43.,  41.,  53.]])
+
+    """
+    m=eye(N, dtype=dtype)
+    if j==None and t==None:
+        raise ValueError("One or more of %s and %s must be set." % \
+            ('j', 't'))
+        return None
+    elif t==None:
+        swap = np.array(m[i])
+        m[i] = m[j]
+        m[j] = swap
+        return m
+    elif j==None:
+        m[i] *= t
+        return m
+    else:
+        m[j] += (t * m[i])
+        return m

-- 
Matt Pagan
matt at pagan.io
PGP: 0xE9284418E360583C

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 801 bytes
Desc: OpenPGP digital signature
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20140324/60e46da2/attachment.sig>