---
title: "High-precision checks for second-tranche SDID weights"
subtitle: "Low-donor and initializer-sensitive problems · fixed default penalty"
format:
  html:
    embed-resources: true
    toc: true
    code-fold: true
    html-math-method: katex
jupyter: python3
---

## Method

Every returned block with at most 20 donors, plus every block whose two QP
initializers differ by more than 10^-5 in L1, is refined at **80-digit precision**.
This covers all six conspicuous initializer disagreements and nearby low-donor
problems. The same outcome values and default noise recipe are used.

A primal active-set algorithm enforces the simplex. It solves the active-face
KKT system directly when the support is small; otherwise a Woodbury identity
moves the solve to the donor dimension. Previously zero periods may re-enter.
Both the QP and default starting vectors are checked. Each solution must satisfy
high-precision feasibility and objective gap below 10^-50, with initializer L1
difference below 10^-45. The remaining blocks retain their double-precision QP
reference; objective and initializer checks alone are not a universal certificate
of weight accuracy under a tiny penalty.

```{python}
#| label: high-precision-solver
import csv, json, math, time
from pathlib import Path
import mpmath as mp
import numpy as np
mp.mp.dps=80
root=Path.cwd().parents[1]
inputs=root/'analyses/sdid-time-weights-tranche2/inputs'
def read(path):
    with open(path) as f:return list(csv.DictReader(f))
manifest={r['fit_id']:r for r in read(inputs/'manifest.csv')}
base=read('fits.csv'); vectors=read('time-weights.csv')
original={r['fit_id']:r for r in base if r['method']=='default'}
initial=read('initialization-sensitivity.csv')
sensitive={r['fit_id'] for r in initial if float(r['l1_difference'])>1e-5}
lookup={}
for r in vectors:lookup.setdefault((r['fit_id'],r['method']),[]).append(float(r['lambda']))

def solve(A,b,start):
    N,n=A.rows,A.cols;rho=mp.mpf('1e-12');tol=mp.mpf('1e-55')
    z=np.asarray(start,dtype=float).copy();z[z<1e-9]=0
    x=mp.matrix([mp.mpf(float(v)) for v in z]);x=x/sum(x)
    S=[i for i in range(n) if x[i]>0];gram={}
    def hij(i,j):
        key=(min(i,j),max(i,j))
        if key not in gram:gram[key]=mp.fsum(A[r,i]*A[r,j] for r in range(N))/N+(rho if i==j else 0)
        return gram[key]
    g=A.T*b/N
    for it in range(1000):
        k=len(S);y=mp.matrix(n,1)
        if k>N:
            # Woodbury solve on the donor dimension, preserving tiny ridge.
            AS=mp.matrix([[A[r,j] for j in S] for r in range(N)])
            K=AS*AS.T+mp.eye(N)*(N*rho)
            v=AS.T*mp.lu_solve(K,b)
            one=mp.matrix([1]*k)
            h=(one-AS.T*mp.lu_solve(K,AS*one))/rho
            mu=(sum(v)-1)/sum(h);ans=v-mu*h
            for i,j in enumerate(S):y[j]=ans[i]
        else:
            K=mp.matrix(k+1);rhs=mp.matrix(k+1,1)
            for i,ii in enumerate(S):
                rhs[i]=g[ii];K[i,k]=1;K[k,i]=1
                for j,jj in enumerate(S):K[i,j]=hij(ii,jj)
            rhs[k]=1;ans=mp.lu_solve(K,rhs)
            for i,ii in enumerate(S):y[ii]=ans[i]
        negative=[i for i in S if y[i]<-tol]
        if negative:
            alpha=min(x[i]/(x[i]-y[i]) for i in negative);x=x+alpha*(y-x)
            for i in S:
                if abs(x[i])<tol:x[i]=0
            S=[i for i in S if x[i]>0];continue
        x=y;gradient=2*(A.T*(A*x-b)/N+rho*x)
        mu=sum(gradient[i] for i in S)/len(S);outside=[i for i in range(n) if i not in S]
        if not outside or min(gradient[i]-mu for i in outside)>=-tol:
            gap=sum(x[i]*gradient[i] for i in range(n))-min(gradient)
            assert min(x)>=-tol and abs(sum(x)-1)<tol and abs(gap)<mp.mpf('1e-50')
            return x,max(mp.mpf(0),gap),it+1
        S.append(min(outside,key=lambda i:gradient[i]))
    raise RuntimeError('Active-set iteration limit')

refined={};receipts=[];export=[]
for r in [r for r in base if r['method']=='qp_uniform' and (int(r['N0'])<=20 or r['fit_id'] in sensitive)]:
    tic=time.time();fid=r['fit_id'];m=manifest[fid];n=int(m['T0'])
    source=np.memmap(inputs/m['matrix_file'],dtype='<f8',mode='r',shape=(int(m['source_N']),int(m['source_T'])))
    selections=read(inputs/m['selection_file'])
    ci=[int(z['matrix_row_1based'])-1 for z in selections if z['role']=='control']
    yf=np.asarray(source[ci,:n+int(m['T1'])])
    Y=[[mp.mpf(float(v)) for v in row] for row in yf]
    N=len(Y);P=len(Y[0])-n
    diffs=[row[t]-row[t-1] for row in Y for t in range(1,n)]
    dm=sum(diffs)/len(diffs);sigma=mp.sqrt(sum((v-dm)**2 for v in diffs)/(len(diffs)-1))
    zeta=mp.mpf('1e-6')*sigma
    assert abs(float(zeta)/(1e-6*float(original[fid]['noise_scale']))-1)<1e-12
    means=[sum(Y[i][t] for i in range(N))/N for t in range(n)]
    target=[sum(row[n:])/P for row in Y];target_mean=sum(target)/N
    A=mp.matrix([[(row[t]-means[t])/sigma for t in range(n)] for row in Y])
    b=mp.matrix([(v-target_mean)/sigma for v in target])
    x,gap,it=solve(A,b,lookup[fid,'qp_uniform'])
    init_diff=mp.mpf(0)
    if True:
        xx,gg,ii=solve(A,b,lookup[fid,'default'])
        init_diff=sum(abs(x[i]-xx[i]) for i in range(n))
        assert init_diff<mp.mpf('1e-45')
    xf=np.asarray([float(v) for v in x]);xf[np.abs(xf)<1e-50]=0;xf/=xf.sum()
    Af=np.asarray(A.tolist(),dtype=float);bf=np.asarray(list(b),dtype=float)
    gradient=2*Af.T@(Af@xf-bf)/N+2e-12*xf
    double_gap=max(0.,float(xf@gradient-gradient.min()))
    default=np.asarray(lookup[fid,'default']);qpu=np.asarray(lookup[fid,'qp_uniform'])
    obj=float(mp.fsum(v*v for v in A*x-b)/N+mp.mpf('1e-12')*mp.fsum(v*v for v in x))
    rr=dict(r);rr.update(method='qp_reference',objective=obj,gap=double_gap,accepted='TRUE',iterations=it,
        iteration_cap='FALSE',solver_status='80-digit active-set KKT verified',raw_violation=0,repair_l1=0,
        seconds=time.time()-tic,l1_from_default=float(np.abs(xf-default).sum()),
        max_change_from_default=float(np.abs(xf-default).max()),
        donor_prediction_rms_change=float(np.sqrt(np.mean((Af@(xf-default))**2))),
        support_1e6=int((xf>1e-6).sum()),effective_periods=float(1/(xf@xf)),
        mass_recent=float(xf[-math.ceil(.2*n):].sum()),mean_position=float(xf@np.linspace(0,1,n)))
    refined[fid]=rr
    template=[w for w in vectors if w['fit_id']==fid and w['method']=='qp_uniform']
    for j,w in enumerate(template):export.append(dict(w,method='qp_reference',lambda_value=float(xf[j]),raw_lambda=float(xf[j])))
    receipts.append({'fit_id':fid,'digits':80,'active_set_steps':it,'mp_gap':mp.nstr(gap,30),
      'initializer_check':True,'mp_initializer_l1':mp.nstr(init_diff,30),
      'l1_from_double_qp':float(np.abs(xf-qpu).sum()),'max_from_double_qp':float(np.abs(xf-qpu).max()),
      'rounded_double_gap':double_gap,'support_1e6':int((xf>1e-6).sum())})
    print(fid,'support',int((xf>1e-6).sum()),'L1 from double QP',float(np.abs(xf-qpu).sum()),flush=True)

def write(name,rows):
    with open(name,'w') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]),lineterminator="\n");w.writeheader();w.writerows(rows)
# Refine low-donor and initializer-sensitive blocks; retain checked QP elsewhere.
# double QP for the remaining papers. The original five variants stay unchanged.
reference=[];reference_weights=[]
for r in [r for r in base if r['method']=='qp_uniform']:
    fid=r['fit_id'];rr=refined.get(fid,dict(r,method='qp_reference'))
    rr['precision_digits']=80 if fid in refined else 16
    reference.append(rr)
for w in vectors:
    if w['method']=='qp_uniform' and w['fit_id'] not in refined:reference_weights.append(dict(w,method='qp_reference'))
for w in export:
    w['lambda']=w.pop('lambda_value');reference_weights.append(w)
write('reference-fits.csv',reference);write('reference-weights.csv',reference_weights)
write('precision-checks.csv',receipts)
Path('precision-verification.json').write_text(json.dumps({'refined_fits':len(refined),'digits':80,
    'initializer_checks':len(refined),'all_mp_gaps_below':'1e-50',
    'all_checked_initializer_l1_below':'1e-45','mpmath_version':mp.__version__},indent=2)+'\n')
```

## Results

```{python}
#| echo: false
import pandas as pd
pd.DataFrame(receipts)
```

The high-precision gap concerns the stated numerical objective, not data accuracy
or statistical uncertainty. The near-zero coefficient makes its minimizer very
sensitive to the precise objective. Large differences from default sparsified
weights can coexist with essentially unchanged fitted donor predictions.

- [High-precision checks](precision-checks.csv)
- [Full optimized reference metadata](reference-fits.csv)
- [Full optimized reference weights](reference-weights.csv)
- [Executable precision-refinement source](precision.qmd)