def leastsq_NxN(x, y, fit_offset=False, perc=None):
    """Solution to least squares: gamma = cov(X,Y) / var(X)
    """
    if perc is not None:
        if not fit_offset and isinstance(perc, (list, tuple)): perc = perc[1]
        weights = csr_matrix(get_weight(x, y, perc)).astype(bool)
        x, y = weights.multiply(x).tocsr(), weights.multiply(y).tocsr()
    else:
        weights = None

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        xx_ = prod_sum_obs(x, x)
        xy_ = prod_sum_obs(x, y)

        if fit_offset:
            n_obs = x.shape[0] if weights is None else sum_obs(weights)
            x_ = sum_obs(x) / n_obs
            y_ = sum_obs(y) / n_obs
            gamma = (xy_ / n_obs - x_ * y_) / (xx_ / n_obs - x_ ** 2)
            offset = y_ - gamma * x_

            # fix negative offsets:
            idx = offset < 0
            gamma[idx] = xy_[idx] / xx_[idx]
            offset = np.clip(offset, 0, None)
        else:
            gamma = xy_ / xx_
            offset = np.zeros(x.shape[1])
    offset[np.isnan(offset)], gamma[np.isnan(gamma)] = 0, 0
    return offset, gamma