<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def _residual_multivariate( self, estimate, unique_R, inds, X, W, template_centers, template_centers_mean_cov, template_widths, template_widths_mean_var_reci, data_sigma):
 """Residual function for estimating centers and widths Parameters estimate : 1D array Initial estimation on centers unique_R : a list of array, Each element contains unique value in one dimension of coordinate matrix R. inds : a list of array, Each element contains the indices to reconstruct one dimension of original cooridnate matrix from the unique array. X : 2D array, with shape [n_voxel, n_tr] fMRI data from one subject. W : 2D array, with shape [K, n_tr] The weight matrix. template_centers: 2D array, with shape [K, n_dim] The template prior on centers template_centers_mean_cov: 2D array, with shape [K, cov_size] The template prior on covariance of centers' mean template_widths: 1D array The template prior on widths template_widths_mean_var_reci: 1D array The reciprocal of template prior on variance of widths' mean data_sigma: float The variance of X. Returns ------- final_err : 1D array The residual function for estimating centers. """


        centers = self.get_centers(estimate)
        widths = self.get_widths(estimate)
        recon = X.size
        other_err = 0 if template_centers is None else (2 * self.K)
        final_err = np.zeros(recon + other_err)
        F = self.get_factors(unique_R, inds, centers, widths)
        sigma = np.zeros((1,))
        sigma[0] = data_sigma
        tfa_extension.recon(final_err[0:recon], X, F, W, sigma)

        if other_err > 0:
            # center error
            for k in np.arange(self.K):
                diff = (centers[k] - template_centers[k])
                cov = from_tri_2_sym(template_centers_mean_cov[k], self.n_dim)
                final_err[recon + k] = math.sqrt(
                    self.sample_scaling *
                    diff.dot(np.linalg.solve(cov, diff.T)))

            # width error
            base = recon + self.K
            dist = template_widths_mean_var_reci *\
                (widths - template_widths) ** 2
            final_err[base:] = np.sqrt(self.sample_scaling * dist).ravel()

        return final_err