Block Vectors and Matrices

Block vectors and matrices (BlockVector and BlockMatrix) provide a mechanism to perform linear algebra operations with very structured matrices and vectors.

When a BlockVector or BlockMatrix is constructed, the number of blocks must be specified.

>>> import numpy as np
>>> from scipy.sparse import coo_matrix
>>> from pyomo.contrib.pynumero.sparse import BlockVector, BlockMatrix
>>> v = BlockVector(3)
>>> m = BlockMatrix(3, 3)

Setting blocks:

>>> v.set_block(0, np.array([-0.67025575, -1.2]))
>>> v.set_block(1, np.array([0.1, 1.14872127]))
>>> v.set_block(2, np.array([1.25]))
>>> v.flatten()
array([-0.67025575, -1.2       ,  0.1       ,  1.14872127,  1.25      ])

The flatten method converts the BlockVector into a NumPy array.

>>> m.set_block(0, 0, coo_matrix(np.array([[1.67025575, 0], [0, 2]])))
>>> m.set_block(0, 1, coo_matrix(np.array([[0, -1.64872127], [0, 1]])))
>>> m.set_block(0, 2, coo_matrix(np.array([[-1.0], [-1]])))
>>> m.set_block(1, 0, coo_matrix(np.array([[0, -1.64872127], [0, 1]])).transpose())
>>> m.set_block(1, 2, coo_matrix(np.array([[-1.0], [0]])))
>>> m.set_block(2, 0, coo_matrix(np.array([[-1.0], [-1]])).transpose())
>>> m.set_block(2, 1, coo_matrix(np.array([[-1.0], [0]])).transpose())
>>> m.tocoo().toarray()
array([[ 1.67025575,  0.        ,  0.        , -1.64872127, -1.        ],
       [ 0.        ,  2.        ,  0.        ,  1.        , -1.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        , -1.        ],
       [-1.64872127,  1.        ,  0.        ,  0.        ,  0.        ],
       [-1.        , -1.        , -1.        ,  0.        ,  0.        ]])

The tocoo method converts the BlockMatrix to a SciPy sparse coo_matrix.

Once the dimensions of a block have been set, they cannot be changed:

>>> v.set_block(0, np.ones(3))
Traceback (most recent call last):
  ...
ValueError: Incompatible dimensions for block 0; got 3; expected 2

Properties:

>>> v.shape
(5,)
>>> v.size
5
>>> v.nblocks
3
>>> v.bshape
(3,)
>>> m.shape
(5, 5)
>>> m.bshape
(3, 3)
>>> m.nnz
12

Much of the BlockVector API matches that of NumPy arrays:

>>> v.sum()
0.62846552
>>> v.max()
1.25
>>> np.abs(v).flatten()
array([0.67025575, 1.2       , 0.1       , 1.14872127, 1.25      ])
>>> (2*v).flatten()
array([-1.3405115 , -2.4       ,  0.2       ,  2.29744254,  2.5       ])
>>> (v + v).flatten()
array([-1.3405115 , -2.4       ,  0.2       ,  2.29744254,  2.5       ])
>>> v.dot(v)
4.781303326558476

Similarly, BlockMatrix behaves very similarly to SciPy sparse matrices:

>>> (2*m).tocoo().toarray()
array([[ 3.3405115 ,  0.        ,  0.        , -3.29744254, -2.        ],
       [ 0.        ,  4.        ,  0.        ,  2.        , -2.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        , -2.        ],
       [-3.29744254,  2.        ,  0.        ,  0.        ,  0.        ],
       [-2.        , -2.        , -2.        ,  0.        ,  0.        ]])
>>> (m - m).tocoo().toarray()
array([[0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.]])
>>> m * v
BlockVector(3,)
>>> (m * v).flatten()
array([-4.26341971, -2.50127873, -1.25      , -0.09493509,  1.77025575])

Accessing blocks

>>> v.get_block(1)
array([0.1       , 1.14872127])
>>> m.get_block(1, 0).toarray()
array([[ 0.        ,  0.        ],
       [-1.64872127,  1.        ]])

Empty blocks in a BlockMatrix return None:

>>> print(m.get_block(1, 1))
None

The dimensions of a blocks in a BlockMatrix can be set without setting a block:

>>> m2 = BlockMatrix(2, 2)
>>> m2.set_row_size(0, 5)
>>> m2.set_block(0, 0, m.get_block(0, 0))
Traceback (most recent call last):
  ...
ValueError: Incompatible row dimensions for row 0; got 2; expected 5.0

Note that operations on BlockVector and BlockMatrix cannot be performed until the dimensions are fully specified:

>>> v2 = BlockVector(3)
>>> v + v2
Traceback (most recent call last):
  ...
NotFullyDefinedBlockVectorError: Operation not allowed with None blocks.
>>> m2 = BlockMatrix(3, 3)
>>> m2 * 2
Traceback (most recent call last):
  ...
NotFullyDefinedBlockMatrixError: Operation not allowed with None rows. Specify at least one block in every row

The has_none property can be used to see if a BlockVector is fully specified. If has_none returns True, then there are None blocks, and the BlockVector is not fully specified.

>>> v.has_none
False
>>> v2.has_none
True

For BlockMatrix, use the has_undefined_row_sizes() and has_undefined_col_sizes() methods:

>>> m.has_undefined_row_sizes()
False
>>> m.has_undefined_col_sizes()
False
>>> m2.has_undefined_row_sizes()
True
>>> m2.has_undefined_col_sizes()
True

To efficiently iterate over non-empty blocks in a BlockMatrix, use the get_block_mask() method, which returns a 2-D array indicating where the non-empty blocks are:

>>> m.get_block_mask(copy=False)
array([[ True,  True,  True],
       [ True, False,  True],
       [ True,  True, False]])
>>> for i, j in zip(*np.nonzero(m.get_block_mask(copy=False))):
...     assert m.get_block(i, j) is not None

Copying data:

>>> v2 = v.copy()
>>> v2.flatten()
array([-0.67025575, -1.2       ,  0.1       ,  1.14872127,  1.25      ])
>>> v2 = v.copy_structure()
>>> v2.block_sizes()
array([2, 2, 1])
>>> v2.copyfrom(v)
>>> v2.flatten()
array([-0.67025575, -1.2       ,  0.1       ,  1.14872127,  1.25      ])
>>> m2 = m.copy()
>>> (m - m2).tocoo().toarray()
array([[0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.]])
>>> m2 = m.copy_structure()
>>> m2.has_undefined_row_sizes()
False
>>> m2.has_undefined_col_sizes()
False
>>> m2.copyfrom(m)
>>> (m - m2).tocoo().toarray()
array([[0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.]])

Nested blocks:

>>> v2 = BlockVector(2)
>>> v2.set_block(0, v)
>>> v2.set_block(1, np.ones(2))
>>> v2.block_sizes()
array([5, 2])
>>> v2.flatten()
array([-0.67025575, -1.2       ,  0.1       ,  1.14872127,  1.25      ,
        1.        ,  1.        ])
>>> v3 = v2.copy_structure()
>>> v3.fill(1)
>>> (v2 + v3).flatten()
array([ 0.32974425, -0.2       ,  1.1       ,  2.14872127,  2.25      ,
        2.        ,  2.        ])
>>> np.abs(v2).flatten()
array([0.67025575, 1.2       , 0.1       , 1.14872127, 1.25      ,
       1.        , 1.        ])
>>> v2.get_block(0)
BlockVector(3,)

Nested BlockMatrix applications work similarly.

For more information, see the API documentation.