code: def image_linear_solve(self, kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps, inv_bool=False):
          
          A = self.linear_response_matrix(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps)
          C_D_response, model_error_list = self.error_response(kwargs_lens, kwargs_ps)
          d = self.data_response
          param, cov_param, wls_model = de_lens.get_param_WLS(A.T, 1 / C_D_response, d, inv_bool=inv_bool)
          kwargs_lens_0 = [kwargs_lens[k] for k in self._idex_lens_list[0]]
          _, _, _, _ = self._imageModel_list[0]._update_linear_kwargs(param, kwargs_lens_0, kwargs_source, kwargs_lens_light, kwargs_ps)
          wls_list = self._array2image_list(wls_model)
          return wls_list, model_error_list, cov_param, param
docs: computes the image (lens and source surface brightness with a given lens model).
          The linear parameters are computed with a weighted linear least square optimization (i.e. flux normalization of the brightness profiles)
          :param kwargs_lens: list of keyword arguments corresponding to the superposition of different lens profiles
          :param kwargs_source: list of keyword arguments corresponding to the superposition of different source light profiles
          :param kwargs_lens_light: list of keyword arguments corresponding to different lens light surface brightness profiles
          :param kwargs_ps: keyword arguments corresponding to "other" parameters, such as external shear and point source image positions
          :param inv_bool: if True, invert the full linear solver Matrix Ax = y for the purpose of the covariance matrix.
          :return: 1d array of surface brightness pixels of the optimal solution of the linear parameters to match the data