def average(A_matrices, B_vectors):
   """ 
   Average over cost function matrices.

   The cost functions contain the three objects of the cost function: A, B, C
   A is a quadratic matrix (97x97), B a vector (d=97), and C is a constant.
   In the end, we are interested in the best-fit charges Q which are the solution to
   Q = A^-1 B

   Arguments:
   A_matrices: a list of NxN matrices
   B_vectors: a list of vectors with len N

   Returns:
   A: the average of all A_matrices.
   B: the average of all B_matrices.
   """

   # Initialize empty
   A = A_matrices[0] * 0
   B = B_vectors[0] * 0

   # Average by adding all objects and dividing by their number
   for index in range(len(A_matrices)):
      A += A_matrices[index]
      B += B_vectors[index]

   # Divide
   number_snapshots = len(A_matrices)
   A /= number_snapshots
   B /= number_snapshots

   return A, B