# Find Cholesky decomposition
cvalue = numpy.linalg.cholesky( covar_matrix )
# Locally store size of matrix
dims = covar_matrix.shape
# Add normal deviate to value, preserving lower triangular
covar_proposed = numpy.multiply( cvalue + random_update, numpy.tri(dims[0]))
# Square and replace
covar_proposed = covar_proposed*transpose(covar_proposed)