code: def estimate_column_scales(
              self,
              X_centered,
              row_scales):
          
          n_rows, n_cols = X_centered.shape
          row_scales = np.asarray(row_scales)
  
          if len(row_scales) != n_rows:
              raise ValueError("Expected length %s, got shape %s" % (
                  n_rows, row_scales.shape,))
  
          column_variances = np.nanmean(
              X_centered ** 2 / (row_scales ** 2).reshape((n_rows, 1)),
              axis=0)
          column_variances[column_variances == 0] = 1.0
          assert len(column_variances) == n_cols, "%d != %d" % (
              len(column_variances),
              n_cols)
          return np.sqrt(column_variances)
docs: column_scale[j] ** 2 =
            mean{i in observed[:, j]}{
              (X[i, j] - row_center[i] - column_center[j]) ** 2
              -------------------------------------------------
                          row_scale[i] ** 2
          }