DCPanelSolver - autorank

[fedemolina]

Hi!

I’ve been experimenting with the spectrum_cut parameter, which is set to 0.002 by default. If I understand correctly, this acts as a form of regularization since the suggested rank is determined by accumulating singular values until they exceed the 1 - spectrum_cut threshold. Essentially, it’s like analyzing the singular value curve, applying a cutoff at a specific threshold, and using that rank to solve the problem.

In other words, we’re defining an energy retention threshold, ensuring that a certain proportion of the total energy is preserved when selecting the rank. That approach seems standard, but I’m curious—why specifically 0.002?

I get that a higher cutoff (e.g., 0.01, keeping 99% of the energy) might discard small but meaningful variations in structured data, though the optimal value could depend on the dataset. On the other hand, a very low cutoff (e.g., 0.0001) might retain too much information, potentially introducing noise or leading to overfitting in treatment periods.

Is there a specific rationale or reference for choosing 0.002? I remember reading that Tibshirani proposed a method for selecting the cutoff, which isn’t the elbow method—does that apply here? (I'm going to look for that paper)

I read in the paper that:
The hyperparameters for all algorithms were tuned using rank 𝑟 ∼ 70 (estimated via cross validation). A test set Ω consisting of 20% of the treated entires is ran- domly sampled and hidden. The test error is then calculated by ‖𝑃Ω(𝑀 + 𝜏𝑍 − 𝑂)‖2 /‖𝑃Ω(𝑂 − 𝑂 ̄)‖2 where 𝑂 ̄ is the mean-value of 𝑂.

But in the code, spectrum_cut is set to a fixed value. Am I missing a theoretical justification for this in the paper?

I looked into it a bit, and it seems like 0.002 might be a commonly used value in the literature. Is that the case?

I experimented with adjusting this parameter because I noticed that the DC method produces strong results, but it also seems to overfit the rows of the outcome matrix.

The first plot displays a denser grid of values, while the last two plots use fewer values.

The final plot specifically shows the mean difference between the observed matrix and the counterfactual matrix. Unlike the previous plots, it does not use multiple Z matrices—just a single one. The calculation follows this approach: tau_vec = np.sum((Y - M) * Z, axis=1) / np.sum(Z, axis=1)

and just to show an individual row of the outcome matrix we can see that by default (maybe not the best colors!) the real curve and the estimated curve are pretty close almost identical sometimes!!

Would it be a good idea to implement an automatic selection method for spectrum_cut? If you agree, I can work on a solution and submit a PR.

[fedemolina]

I took another look at the current implementation, and I think there are a few separate questions mixed together here:

1. What does spectrum_cut=0.002 actually do?
2. Is the current auto-rank formula implemented correctly?
3. Is an explained-energy cutoff the right way to choose the rank for DC-PR?
4. Would cross-validation be a better long-term approach?

Current implementation

Right now, DCPanelSolver.fit(...) defaults to:

def fit(self, suggest_r=None, auto_rank=True, spectrum_cut=0.002, method='auto', method_non_neg=None):
    if suggest_r is not None:
        M, tau, std = self.DC_PR_with_suggested_rank(...)
    elif auto_rank:
        M, tau, std = self.DC_PR_auto_rank(
            spectrum_cut=spectrum_cut,
            method=method,
            method_non_neg=method_non_neg,
        )

and DC_PR_auto_rank does:

def DC_PR_auto_rank(self, spectrum_cut=0.002, method='convex', method_non_neg=None):
    s = np.linalg.svd(self.O, full_matrices=False, compute_uv=False)
    suggest_r = np.sum(np.cumsum(s**2) / np.sum(s**2) <= 1-spectrum_cut)
    return self.DC_PR_with_suggested_rank(
        suggest_r=suggest_r,
        method=method,
        method_non_neg=method_non_neg,
    )

So spectrum_cut is currently an explained-energy cutoff over squared singular values, not raw singular values. In other words, it keeps components until the cumulative Frobenius energy gets near 1 - spectrum_cut.

For the default value:

spectrum_cut = 0.002

the intended interpretation seems to be something like: choose a rank that explains roughly 99.8% of the observed matrix energy.

Possible off-by-one issue

I think the current formula may be off by one relative to the usual interpretation of an energy threshold.

The current code is:

suggest_r = np.sum(np.cumsum(s**2) / np.sum(s**2) <= 1-spectrum_cut)

This counts how many cumulative energy values are still below or equal to the threshold. But if the goal is to select the first rank that reaches the threshold, the usual implementation would be closer to:

energy = np.cumsum(s**2) / np.sum(s**2)
suggest_r = np.searchsorted(energy, 1 - spectrum_cut, side="left") + 1

with a clamp to ensure:

1 <= suggest_r <= len(s)

The current implementation can even return 0 for reasonable spectrum_cut values if the first singular component already explains more than the threshold.

For example, in the California smoking example, the first singular component explains about 99.25% of the squared singular-value energy. Then:

spectrum_cut = 0.01  # threshold = 99%

would produce:

current formula: rank 0
first rank reaching threshold: rank 1

That seems undesirable.

With the current default:

spectrum_cut = 0.002  # threshold = 99.8%

the current formula also seems to select the rank just before the threshold is crossed, rather than the first rank that reaches the threshold.

So independent of the broader statistical question, I think there may be a small correctness fix here.

Why 0.002?

I do not see an obvious theoretical justification in the code or docs for this exact value.

The docs currently say that DC_PR_auto_rank:

helps to find the proper rank for M (but not very stable, and may be updated later)

That matches what I am seeing: 0.002 behaves like a heuristic. It may work reasonably in some examples, but it is not obviously tied to a theorem in the DC-PR paper or to the cross-validation procedure described there.

So I think we should be careful documenting it as a heuristic default rather than a theoretically tuned parameter.

Rank selection on raw O vs adjusted outcome

Another thing worth discussing: DC_PR_auto_rank currently computes the singular values of self.O directly:

s = np.linalg.svd(self.O, full_matrices=False, compute_uv=False)

But the model is closer to:

O ≈ M + tau Z

So the treatment signal can affect the spectrum of O. This may matter when treatment effects are large or treatment patterns are structured.

One possible improvement would be to estimate a preliminary treatment effect first:

pre_tau = self.solve_tau(self.O)
O_adj = self.O - np.tensordot(self.Z, pre_tau, axes=([2], [0]))

and then choose the rank based on O_adj rather than raw O.

That would be more consistent with the logic already used inside DC_PR_with_suggested_rank, where pre_tau is computed before the convex path chooses the regularization trajectory.

This is still heuristic, but it may be a better heuristic than selecting rank directly from raw outcomes.

Cross-validation

The paper’s approach seems closer to cross-validating the rank or regularization parameter by hiding entries and evaluating reconstruction/prediction error.

That would probably be better than a fixed spectrum_cut, but I think it is a larger design change.

The main complication is that the current DC solver does not appear to support a mask or weighted loss. If we simply fit on the full matrix and then evaluate “held-out” entries, that would leak information because those entries were used in fitting.

A proper CV version would likely require one of these:

• add an observation mask / sample-weight matrix to the DC objective;
• implement a masked version of the alternating solver;
• or design a validation scheme over units/time blocks that avoids leakage.

So I would probably not start with full CV unless we are ready to modify the solver more deeply.

Suggested incremental path

I think a good first PR could be small and focused:

1. Add a helper function for rank selection from singular values, e.g.

def rank_from_spectrum(s, spectrum_cut):
    ...

2. Validate spectrum_cut, e.g.

0 < spectrum_cut < 1

3. Use np.searchsorted to choose the first rank that reaches the target energy.

4. Clamp the rank to [1, len(s)].

5. Add unit tests with synthetic singular value arrays to verify:

   • no rank 0;
   • exact threshold behavior;
   • small spectrum_cut gives larger ranks;
   • large spectrum_cut gives smaller ranks.

6. Document that this is an explained-energy heuristic.

A possible implementation sketch:

def rank_from_spectrum(s, spectrum_cut=0.002):
    if not 0 < spectrum_cut < 1:
        raise ValueError("spectrum_cut must be between 0 and 1.")

    s = np.asarray(s)
    energy_total = np.sum(s**2)
    if energy_total <= 0:
        return 1

    cumulative_energy = np.cumsum(s**2) / energy_total
    target = 1 - spectrum_cut
    rank = np.searchsorted(cumulative_energy, target, side="left") + 1
    return int(np.clip(rank, 1, len(s)))

Then:

def DC_PR_auto_rank(self, spectrum_cut=0.002, method='convex', method_non_neg=None):
    s = np.linalg.svd(self.O, full_matrices=False, compute_uv=False)
    suggest_r = rank_from_spectrum(s, spectrum_cut=spectrum_cut)
    return self.DC_PR_with_suggested_rank(
        suggest_r=suggest_r,
        method=method,
        method_non_neg=method_non_neg,
    )

A slightly more opinionated version could select rank from the preliminary treatment-adjusted matrix:

pre_tau = self.solve_tau(self.O)
O_adj = self.O - np.tensordot(self.Z, pre_tau, axes=([2], [0]))
s = np.linalg.svd(O_adj, full_matrices=False, compute_uv=False)
suggest_r = rank_from_spectrum(s, spectrum_cut=spectrum_cut)

But I think that should probably be discussed separately because it changes existing behavior more than just fixing the threshold mapping.

Open design question

Maybe we should distinguish two modes explicitly:

solver.fit(auto_rank=True, rank_selection="energy")
solver.fit(auto_rank=True, rank_selection="cv")

where "energy" keeps the current lightweight heuristic, and "cv" can later implement the paper-style validation procedure once the solver supports masked fitting.

For now, I think the safest first step is to fix/document the energy heuristic and add tests. Then we can revisit cross-validation as a larger follow-up.