Let us start with the matrices stored with the ndarray class. We accomplish scalar multiplication with the * operator, and the matrix addition with the + operator. But for matrix multiplication we will need the instance method dot() or the numpy.dot function, since the operator * is reserved for element-wise multiplication:

In [54]: 2*A
Out[54]:
array([[ 2,  4],
       [ 8, 32]])
In [55]: A + 2*A
Out[55]:
array([[ 3,  6],
       [12, 48]])
In [56]: A.dot(2*A)       In [56]: np.dot(A, 2*A)
Out[56]:                  Out[56]:
array([[ 18,  68],        array([[ 18,  68],
       [136, 528]])              [136, 528]])
In [57]: A.dot(B)
ValueError: objects are not aligned
In [58]: B.dot(A)         In [58]: np.dot(B, A)
Out[58]:                  Out[58]:
array([[ -9, -34],        array([[ -9, -34],
       [  0,   0],               [  0,   0],
       [  9,  34]])              [  9,  34]])

The matrix class makes matrix multiplication more intuitive: the operator * can be used instead of the dot() method. Note also how matrix multiplication between different instance classes ndarray and a matrix is always casted to a matrix instance class:

In [59]: C * B
ValueError: shapes (2,2) and (3,2) not aligned: 2 (dim 1) != 3 (dim 0)
In [60]: B * C
Out[60]: 
matrix([[ -9, -34],
        [  0,   0],
        [  9,  34]])

For sparse matrices, both scalar multiplication and addition work well with the obvious operators, even if the two sparse classes are not the same. Note the resulting class casting after each operation:

In [61]: S_10_coo = spsp.rand(5, 5, density=0.1, format='coo')
In [62]: S_25_lil + S_10_coo
Out[62]: <5x5 sparse matrix of type '<type 'numpy.float64'>'
        with 8 stored elements in Compressed Sparse Row format>
In [63]: S_25_lil * S_10_coo
Out[63]: <5x5 sparse matrix of type '<type 'numpy.float64'>'
        with 4 stored elements in Compressed Sparse Row format>

In [64]: S_100_coo = spsp.rand(2, 2, density=1, format='coo')
In [65]: np.dot(A, S_100_coo)
Out[66]:
array([[ <2x2 sparse matrix of type '<type 'numpy.float64'>'
  with 4 stored elements in COOrdinate format>,
        <2x2 sparse matrix of type '<type 'numpy.float64'>'
  with 4 stored elements in COOrdinate format>],
       [ <2x2 sparse matrix of type '<type 'numpy.float64'>'
  with 4 stored elements in COOrdinate format>,
        <2x2 sparse matrix of type '<type 'numpy.float64'>'
  with 4 stored elements in COOrdinate format>]], dtype=object)
In [67]: A * S_100_coo
Out[68]:
array([[  1.81 ,   1.555],
       [ 11.438,  11.105]])

The traces of a matrix are the sums of the elements on the diagonals (assuming always increasing indices in both dimensions). For generic matrices, we compute them with the instance method trace(), or with the function numpy.trace:

In [69]: A.trace()        In [71]: C.trace()
Out[69]: 17               Out[71]: matrix([[17]])
In [70]: B.trace()        In [72]: np.trace(B, offset=-1)
Out[70]: -1               Out[72]: 2

In order to compute the determinant of generic square matrices, we need the function det in the module scipy.linalg:

In [73]: spla.det(C)
Out[73]: 8.0

Transposes can be computed with any of the two instance methods transpose() or T, for any of the two classes of generic matrices:

In [74]: B.transpose()        In [75]: C.T
Out[74]:                      Out[75]: 
array([[-1,  0,  1],          matrix([[ 1,  4],
       [-2,  0,  2]])                 [ 2, 16]])

Hermitian transpose can be computed for the matrix class with the instance method H:

In [76]: D = C * np.diag((1j,4)); print D     In [77]: print D.H
[[  0.+1.j   8.+0.j]                          [[  0.-1.j   0.-4.j]
 [  0.+4.j  64.+0.j]]                          [  8.-0.j  64.-0.j]]

Inverses of non-singular square matrices are computed for the ndarray class with the function inv in the module scipy.linalg. For the matrix class, we may also use the instance method I. For non-singular square sparse matrices, we may use the function inv in the module scipy.sparse.linalg.

In [78]: E = spsp.rand(512, 512, density=1).todense()
In [79]: S_100_csc = spsp.rand(512, 512, density=1, format='csc')
In [80]: %timeit E.I
10 loops, best of 3: 28.7 ms per loop
In [81]: %timeit spspla.inv(S_100_csc)
1 loops, best of 3: 1.99 s per loop

/scipy/sparse/linalg/dsolve/linsolve.py:88: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  warn('spsolve requires A be CSC or CSR matrix format', SparseEfficiencyWarning)
/scipy/sparse/linalg/dsolve/linsolve.py:103: SparseEfficiencyWarning: solve requires b be CSC or CSR matrix format

The Moore-Penrose pseudo-inverse can be computed for any kind of matrix (not necessarily square) with either routines the pinv or the pinv2 in the module scipy.linalg. The first method, pinv, resorts to solving a least squares problem to compute the pseudo-inverse. The function pinv2 computes the pseudo-inverse by a method based on singular value decompositions. For Hermitian matrices, or matrices that are symmetric with no complex coefficients, we also have a third function called pinvh, which is based on eigenvalue decompositions.

It is known that in the case of square non-singular matrices, the inverse and pseudo-inverse are the same. This simple example shows the times of computation of the inverses of a large generic symmetric matrix with the five methods described:

In [82]: F = E + E.T     # F is symmetric
In [83]: %timeit F.I
1 loops, best of 3: 24 ms per loop
In [84]: %timeit spla.inv(F)
10 loops, best of 3: 28 ms per loop
In [85]: %timeit spla.pinvh(E)
1 loops, best of 3: 120 ms per loop
In [86]: %timeit spla.pinv2(E)
1 loops, best of 3: 252 ms per loop
In [87]: %timeit spla.pinv(F)
1 loops, best of 3: 2.21 s per loop

For generic matrices, we have seven different standard norms in scipy.linalg. We can summarize them in the following table:

In [88]: [spla.norm(A,s) for s in (np.inf,-np.inf,-1,1,-2,2,'fro')]
Out[88]: [20, 3, 5, 18, 0.48087417361008861, 16.636368595013604, 16.643316977093239]

In [89]: spla.norm(S_100_csc.todense(), 1) - \
   ....: spspla.onenormest(S_100_csc)
Out[89]: 0.0

Note the subtle difference between the norm computations from SciPy and NumPy. For example, in the case of the Frobenius norm, scipy.linalg.norm is based directly on the BLAS function called NRM2, while numpy.linalg.norm is equivalent to a purely straightforward computation of the form sqrt(add.reduce((x.conj() * x).real)). The advantage of the code based on BLAS, besides being much faster, is clear when some of the data is too large or too small in single-precision arithmetic. This is shown in the following example:

In [89]: a = np.float64([1e20]); \
   ....: b = np.float32([1e20])
In [90]: [np.linalg.norm(a), spla.norm(a)]
Out[90]: [1e+20, 1e+20]
In [91]: np.linalg.norm(b)
[...]/numpy/linalg/linalg.py:2056: RuntimeWarning: overflow encountered in multiply
  return sqrt(add.reduce((x.conj() * x).real, axis=None))
Out[91]: inf
In [92]: spla.norm(b)
Out[92]: 1.0000000200408773e+20

This brings us inevitably to a discussion about the computation of the condition number of a non-singular square matrix A. This value measures how much the output of the solution to the linear equation A * x = b will change when we make small changes to the input argument b. If this value is close to one, we can rest assured that the solution is going to change very little (we say then that the system is well-conditioned). If the condition number is large, we know that there might be issues with the computed solutions of the system (and we say then that it is ill-conditioned).

The computation of this condition number is performed by multiplying the norm of A with the norm of its inverse. Note that there are different condition numbers, depending on the norm that we choose for the computation. These values can also be computed for each of the pre-defined norms with the function numpy.linalg.cond, although we need to be aware of its obvious limitations.

In [93]: np.linalg.cond(C, -np.inf)
Out[93]: 1.875