\section{Introduction}
\paragraph{Motivation.}
Least mean squares solvers are fundamental tools in all the data science fields such as machine learning, computer science and statistics. They are also the building blocks of more involved techniques such as deep learning and signal processing~\cite{mandic2004generalized,widrow1977stationary}.  As explained in~\cite{maalouf2019fast}, this family include Singular Value Decomposition (SVD), Principle Component Analysis (PCA), linear regression, Lasso and Ridge regression, Elastic net, and many more~\cite{golub1971singular,jolliffe2011principal,hoerl1970ridge,seber2012linear,zou2005regularization,tibshirani1996regression,safavian1991survey}. First closed form solutions for problems such as linear regression were published by e.g. Pearson~\cite{pearson1900x} around 1900 but were probably known before. Nevertheless, today they are still used extensively as building blocks in both academy and industry for normalization~\cite{liang2013distributed,kang2011scalable,afrabandpey2016regression}, spectral clustering~\cite{peng2015robust}, graph theory~\cite{zhang2018understanding}, prediction~\cite{copas1983regression,porco2015low}, dimensionality reduction~\cite{laparra2015dimensionality}, feature selection~\cite{gallagher2017cross} and many more; see more examples in~\cite{golub2012matrix}.

Important special case is the low rank approximation of an $n\times d$ real matrix $A$ that can be computed via $k$-SVD (Singular Value Decomposition). Which is the linear (non-affine) $k$-dimensional subspace that minimizes its sum of squared distances over the rows of $A$ for a given integer $k\geq1$, i.e.,
\[
\argmin_{X\in\ensuremath{\mathbb{R}}^{d\times k},X^TX=I } \sum_{i=1}^n \norm{A_i - A_iXX^T}^2_2,
\]
where $A_i$ is the $i$th row of the matrix $A$ for every integer $i\in \br{1,\cdots , n}$.

More generally, the $k$-PCA is the \emph{affine} $k$-subspace that minimizes the sum of squared distances from the rows of $A$ to it over every $k$-subspace that may be translated from the origin of $R^d$. Formally, an affine $k$-subspace is represented by an orthogonal matrix $X\in \ensuremath{\mathbb{R}}^{d\times k}$ and a vector $\ell\in \ensuremath{\mathbb{R}}^d$ that represent the translation of the subspace from the origin. Hence, we wish to compute:
\[
\argmin_{\substack{ \ell\in \ensuremath{\mathbb{R}}^d \\ X\in\ensuremath{\mathbb{R}}^{d\times k} , X^TX=I }} \sum_{i=1}^n \norm{(A_i-\ell^T) - (A_i-\ell^T)XX^T}^2_2.
\]

Finally we have the Least-Mean-Squares solvers that gets as input an $n\times d$ real matrix $A$, and another $n$-dimensional real vector $b$ (possibly the zero vector), and aims to minimize the sum of squared distances from the rows (points) of $A$ to some hyperplane that is represented by its normal or vector of $d$ coefficients $x$, that is constrained to be in a given set $X\subseteq\ensuremath{\mathbb{R}}^d$:
\[
\argmin_{x\in X}f(\norm{Ax-b}_2)+g(x).
\]
Here, $g$ is called a \emph{regularization term}. For example: in linear regression $X=\ensuremath{\mathbb{R}}^d$, $f(x)=x^2$ and $g(x)=0$ for every $x\in X$. In Lasso $f(y) = y^2$ and $g(y) = \alpha \cdot \norm{x}_1$ for every $y\in \ensuremath{\mathbb{R}}^d$ and $\alpha>0$.
Such LMS solvers can be computed via the covariance matrix $A^TA$. For example, the solution to linear regression of minimizing $\norm{Ax-b}_2$ is $(A^TA)^{-1}A^Tb$.
\subsection{Coresets}
 For a huge amount of data, those algorithms/sovlers are much time consuming: while in theory the running time is usually $O(nd^2)$, the constants that are hidden in the $O(\cdot)$ notation are significantly large.
 Another problem with such algorithms/sovlers is that we may not be able to use them for big data  on standard machines, since there is no enough memory to provide the relevant computations.

A modern tool to handle these type of problems, is a data summarization for the input that is sometimes called \emph{coresets}. Coresets also allow us to boost the running time of those algorithms/solvers while using less memory.

As explained at~\cite{unpublishedDanny}, coresets are especially useful to (a) learn unbounded streaming data that cannot fit into main memory, (b) run in parallel on distributed data among thousands of machines, (c) use low communication between the machines, (d) apply real-time computations on a device, (e) handle privacy and security issues, (f) compute constrained optimization on a coreset that was constructed independently of these constraints and of course boost there running time.

In the context of the $k$-SVD problem, an $\varepsilon$-coreset for a matrix $A\in\ensuremath{\mathbb{R}}^{n\times d}$ is a matrix $C\in \ensuremath{\mathbb{R}}^{m\times d}$ where $m \ll n$, which guarantees that the sum of the squared distances from any linear (non-affine) $k$-subspace to the rows of $C$ will be approximately equal to the  sum of the squared distances from the same $k$-subspace to the rows of $A$, up to a $(1\pm\varepsilon)$ multiplicative factor, i.e., for any matrix $X\in\ensuremath{\mathbb{R}}^{d\times k}$  such that $X^TX=I$ we have,
$$ \abs{\sum_{i=1}^n \norm{A_i - A_iXX^T}^2_2 -  \sum_{i=1}^m \norm{C_i - C_iXX^T}^2_2} \leq \varepsilon \sum_{i=1}^n \norm{A_i - A_iXX^T}^2_2.$$

In the $k$-PCA problem, an $\varepsilon$-coreset for the matrix $A$ is a matrix $C\in \ensuremath{\mathbb{R}}^{m\times d}$ such that for every vector $\ell \in \ensuremath{\mathbb{R}}^d$  and a matrix $X\in\ensuremath{\mathbb{R}}^{d\times k}$ where $X^TX=I$ we have:
\[
\abs{ \sum_{i=1}^n \norm{(A_i-\ell^T) - (A_i-\ell^T)XX^T}^2_2 - \sum_{i=1}^m\norm{(C_i-\ell^T) - (C_i-\ell^T)XX^T}^2_2} \leq \varepsilon\sum_{i=1}^n \norm{(A_i-\ell^T) - (A_i-\ell^T)XX^T}^2_2.
\]

The dimension of the subspace $k$ is a crucial parameter, and of course may determines the size of the coreset and the time complexity of the algorithm.
Such coresets are useful, for example, for many NLP applications in which a Word Embedding model needed to be produced out of a large database, see \cite{mikolov2013distributed,mikolov2013efficient,pennington2014glove}.

Considering  the least means squared problems. Given a matrix  $A\in \ensuremath{\mathbb{R}}^{n\times d}$ and a $n$-dimensional real vector $b$, a coreset for the pair $(A,b)$ in the LMS problem that is represented by the functions $f$ and $g$ as explained before, is a matrix $C\in \ensuremath{\mathbb{R}}^{m\times d}$ where $m \ll n$ and a vector $y\in \ensuremath{\mathbb{R}}^{m}$, such that  for every hyperplane that is represented by its normal or vector of $d$ coefficients $x$ we have,
$$\abs{f(\norm{Ax-b}_2)- f(\norm{Cx-y}_2)} \leq \varepsilon \big( f(\norm{Ax-b}_2)+g(x) \big).$$
Usually $g$ is  a non-negative function (e.g., $g(x) = \norm{x}_2^2$ in Ridge, $g(x) = \norm{x}_1$ in Lasso, and $g(x)=0$ in linear regression). Hence for those cases, $(C,y)$ is a coreset for $(A,b)$  if it satisfies:
\[
\abs{f(\norm{Ax-b}_2)- f(\norm{Cx-y}_2) }\leq \varepsilon f(\norm{Ax-b}_2),
\]
for every $x\in \ensuremath{\mathbb{R}}^d$.

\subsection{Coreset constructions}
One type of coresets, sometimes called \emph{sketch}, consists on linear combinations of the input points. These coresets use techniques such as Random projections~\cite{DBLP:journals/corr/CohenEMMP14}, JL-Lemma~\cite{sarlos2006improved},SVD~\cite{DBLP:journals/corr/abs-1807-04518} etc.
However, in this paper we consider only coresets that are \emph{subset} of their input points (rows of the input matrix), up to a multiplicative weight (scaling).

As explained in~\cite{unpublishedDanny} and~\cite{DBLP:journals/corr/BargerF15}, the advantages of such coresets are: (i) preserved sparsity of the input, (ii) interpretability, (iii) coreset may be used (heuristically) for other problems, (iv)~less numerical issues that occur when non-exact linear combination of points are used.

\paragraph{Sensitivity sampling.}
Over the recent decades many algorithms were suggested to compute such coresets.
One of the common technique, both in theory and practice, that yields fast and provably good coresets is the approach of non-uniform sampling, sensitivity sampling~\cite{langberg2010universal,DBLP:journals/corr/BravermanFL16}. The sensitivity of a row $p$ in the input matrix $A$ is a number $s(p)\in [0,1]$ that represents how much this row is `important' in this dataset with respect to the desired optimization problem. The motivation for defining sensitivity is the following. Suppose that we have an upper bound $s'(p)\geq s(p)$ for every row $p$ in the matrix $A$, and we use it to sample $m\geq 1$ rows from $A$ where every row $p$ is picked with probability $\mathrm{pr}(p)$ that is proportional to $s'(p)$ and is assigned a weight of $w(p)=\frac{1}{m\cdot pr(p)}$. Then sampling such $\log(d)T/\varepsilon^2$ i.i.d. rows would yield an $\varepsilon$-coreset where $T=\sum_{i=1}^n s'(A_i)$ is called the total sensitivity bound ($A_i$ is the $i$th row of $A$ for every integer $i\in \br{1,\cdots,n}$ as previously defined).

\paragraph{Sensitivity bounds. }One of the main challenges in constructing coresets is to bound the corresponding sensitivity of each point, which is $s(\cdot)$ in~\eqref{sensreduction}.
While $1$ is a trivial bound for $s(\cdot)$, it would give a coreset of size $|C|=O(n)$ as $T=n$ in this case. In the $k$-SVD problem, for the case of $k=d-1$, as shown at~\cite{yang2017weighted} the sensitivity is also known as leverage score and can be easily bounded by $s'(p)=\norm{u}_2^2$ where $u$ is the corresponding row for $p$ in the matrix $U$ such that $A=UDV^T$ is the thin SVD of the matrix $A$; see Definition~\ref{thinsvddef}, and the sum of sensitivities is exactly $T$. It is easy to prove that this bound is tight in the sense that $s'(p)=s(p)$. For the case $k\leq d-2$, sensitivity bounds are also known whose total sensitivity is $T=O(d)$, by projecting the points on an optimal (or approximated) $k$-subspace an computing the sensitivity of the projected point as shown at~\cite{DBLP:journals/corr/abs-1209-4893}. However, unlike the previous case, these bounds are not tight, as proved in the experimental results of this paper. In the $k$-PCA problem, a tight sensitivity bound for the case $k=0$ (i.e., the $1$-mean problem) was suggested at~\cite{tremblay2018determinantal}, however there is no tight bound for the other cases.

\subsection{Our contribution}
In this work we suggest:
\begin{enumerate} [(i)]
\setlength{\itemsep}{-1pt}
\item  the first algorithm that computes tight sensitivity bounds for the family of (non-affine) $k$-subspaces; see Algorithm~\ref{alg:tight1}. The algorithm is iterative and returns the exact sensitivity $s(p)$ up to arbitrarily small constant. The convergence rate is linear.
\item generalization of the above algorithm for the family of affine $k$-dimensional subspaces of $\ensuremath{\mathbb{R}}^d$. This is by reduction to a problem of computing sensitivity bounds of a new set of points in $\ensuremath{\mathbb{R}}^{d+1}$ for the family of (non-affine) $(k+1)$- subspaces in $\ensuremath{\mathbb{R}}^{d+1}$. See Theorem~\ref{lemma:pca}.
\item experimental results on real-world datasets, including the English Wikipedia documents-term matrix, that show that our bounds provide significantly smaller and data-dependent coresets also in practice.
\item full open source code.
\end{enumerate}
While our sensitivity bounds are tight for the family of affine (or non-affine) $k$-subspaces in $\ensuremath{\mathbb{R}}^d$ they are no longer tight if we consider only subset of this family of subspaces. Nevertheless, they provide better upper bounds for these problems or query spaces, compared to existing upper bounds that also ignore these constraints and regularization terms. More precisely, the worst case sensitivity is $O(k)$ in both cases, but our bounds are tighter.

\section{Preliminaries}
In the this section we give our notations and definitions that will be used through the paper. We also explain the relation between the notion of total sensitivity and coreset size while relying on Theorem~5.5 in~\cite{DBLP:journals/corr/BravermanFL16}.

\begin{paragraph}{Notations.}For integers $d,n\geq 1$, we denote by $0_d$ the origin of $\ensuremath{\mathbb{R}}^d$. The set $\ensuremath{\mathbb{R}}^{n\times d}$ denote the union over every $n\times d$ real matrix. For a matrix $A\in\ensuremath{\mathbb{R}}^{n\times d}$ the Frobenius norm $\norm{A}_F$ is the squared root of its sum of squared entries, and $Tr(A)$ denotes its trace. A \emph{weighted set} of $n$ points in $\ensuremath{\mathbb{R}}^d$ is a pair $(P,w)$ where $P=\br{p_1,\cdots,p_n}$ is an ordered set in $\ensuremath{\mathbb{R}}^d$, and $w: P\to [0,\infty)$ is called a \emph{weight function}.

For an integer $k \in \br{0,\cdots,d-1}$, a \emph{$k$-subspace} is a shorthand for a $k$-dimensional linear (non-affine) subspace of $\ensuremath{\mathbb{R}}^d$ (i.e., it contains the origin). An \emph{affine $k$-subspace} ($k$-flat) is a translation of a $k$-subspace, i.e.,  that may not contain the origin.
For every point $p\in \ensuremath{\mathbb{R}}^d$ and an affine $k$-subspace $S$ of $\ensuremath{\mathbb{R}}^d$, we define $\mathrm{proj}(p,S)=\argmin_{x\in S} \norm{p-x}_2$ to be the projection of the point $p$ onto the affine $k$-subspace $S$ and $\mathrm{dist}(p,S)=\min_{x\in S}\norm{p-x}_2=\norm{p-\mathrm{proj}(p,S)}_2$ to be the Euclidean distance between the point $p$ to its closest point on $S$. This distance to the power of $z\geq1$ is denoted by $D_z(p,S) =\mathrm{dist}^z(p,S)$ and for brevity, we define $D(p,S) =D_2(p,S) =\mathrm{dist}^2(p,S)$.
\end{paragraph}

\begin{definition} [\emph{Additive $\eps$-approximation}]
Let $d$ be an integer, $\varepsilon\in (0,1)$ be an error parameter, $f:\ensuremath{\mathbb{R}}^d \to \ensuremath{\mathbb{R}}$ be a function and $s\in \ensuremath{\mathbb{R}}$ be a real number. We call $s$ an \emph{additive $\eps$-approximation} for $f$ if and only if
$$\sup_{x\in \ensuremath{\mathbb{R}}^d}f(x) \leq s \leq \sup_{x\in \ensuremath{\mathbb{R}}^d}f(x)+\varepsilon .$$
\end{definition}

\begin{definition} [\emph{Thin SVD}] \label{thinsvddef}
Let $n,d$ be two integers. Let $\mathbf{A}\in \ensuremath{\mathbb{R}}^{n\times d}$ be a matrix and let the integer $r$ be its rank. We call $\mathbf{A}=UDV^T$ the thin Singular Value Decomposition of $\mathbf{A}$. That is, $U\in \ensuremath{\mathbb{R}}^{n\times r}$, $V\in \ensuremath{\mathbb{R}}^{d\times r}$, $U^TU=I$, $V^TV=I$ and $D$ is a diagonal matrix.
\end{definition}

\begin{definition} [\textbf{Definition 4.2 in~\cite{DBLP:journals/corr/BravermanFL16}}] \label{def:querySpace}
Let $(P,w)$ be a weighted set of $n$ points in $\ensuremath{\mathbb{R}}^d$. Let $Q$ be a function that maps every set $C\subseteq P$ to a corresponding set $Q(C)$, such that $Q(T) \subseteq Q(C)$ for every $T\subseteq C$. Let $f:P\times Q(P) \to \ensuremath{\mathbb{R}}$ be a cost function. The tuple $(P,w,Q,f)$ is called a \emph{query space}.
\end{definition}

\begin{definition} [\textbf{Definition 4.5 in~\cite{DBLP:journals/corr/BravermanFL16}}] \label{def:VC}
For a query space $(P,w,Q,f)$, $q\in Q(P)$ and $r\in [0,\infty)$ we define
\[
\mathrm{range}(q,r) = \br{p\in P \mid w(p)\cdot f(p,q) \leq r}.
\]
The dimension of $(P,w,Q,f)$ is the smallest integer $d'$ such that for every $C\subseteq P$ we have
\[
\left| \br{\mathrm{range}(q,r) \mid q\in Q(C), r\in [0,\infty)} \right| \leq |C|^{d'}.
\]
\end{definition}

\begin{theorem} [\textbf{Theorem 5.5 in~\cite{DBLP:journals/corr/BravermanFL16}}]\label{braverman}
Let $(P,w,Q,f)$ be a query space; see Definition~\ref{def:querySpace}, where $f$ is a non-negative function. Let $s:P\to [0,\infty)$ such that
\[
\sup_{q\in Q} \frac{w(p)f(p,q)}{\sum_{p\in P} w(p)f(p,q)} \leq s(p),
\]
for every $p\in P$ and $q\in Q(P)$ such that the denominator is non-zero. Let $t = \sum_{p\in P} s(p)$ and let $d'$ be the dimension of the query space $(P,w,Q,f)$; See Definition~\ref{def:VC}. Let $c \geq 1$ be a sufficiently large constant and let $\varepsilon, \delta \in (0,1)$. Let $C$ be a random sample of
\[
|C| \geq \frac{ct}{\varepsilon^2}\left(d'\log{t}+\log{\frac{1}{\delta}}\right)
\]
points from $P$, such that $p$ is sampled with probability $s(p)/t$ for every $p\in P$. Let $u(p) = \frac{t\cdot w(p)}{s(p)|C|}$ for every $p\in C$. Then, with probability at least $1-\delta$, for every $q\in Q$ it holds that
\[
(1-\varepsilon)\sum_{p\in P} w(p)\cdot f(p,q) \leq \sum_{p\in C} u(p)\cdot f(p,q) \leq (1+\varepsilon)\sum_{p\in P} w(p)\cdot f(p,q).
\]
\end{theorem}
\begin{paragraph} {Smaller total sensitivity implies smaller coreset size.}
Given $(P,w)$ a weighted set of $n$ points in $\ensuremath{\mathbb{R}}^d$ and given also the sensitivity $s:P\to [0,\infty)$ of each point as defined in~\ref{braverman}, in order to obtain a coreset that guarantees $(1\pm\varepsilon)$ multiplicative error with probability at least $1-\delta$ we have to sample $O(\frac{ct}{\varepsilon^2}\left(d'\log{t}+\log{\frac{1}{\delta}}\right))$ points from $P$ where $t=\sum_{i=1}^n s(p)$ is the sum of sensitivity over all the points in $P$. Thus the smaller the sensitivity bound $s(p)$ of each point $p\in P$ the smaller the total sensitivity $t$ and the smaller is the size of the coreset needed.

\end{paragraph}
\section{Sensitivity of  Non-affine $k$-subspaces}
Let $S'$ be a (non-affine) $k$-subspace of $\ensuremath{\mathbb{R}}^d$. Every such subspace $S'$ corresponds to a column space of a matrix $X\in \ensuremath{\mathbb{R}}^{d\times k}$ whose columns are orthonormal ($X^TX=I$). Let $(P,w)$ be a weighted set of $n\geq 1$ points in $\ensuremath{\mathbb{R}}^d$.
Let $\mathbf{P}\in \ensuremath{\mathbb{R}}^{n\times d}$ denote the matrix whose $i$th row is the $i$th point of $P$ multiplied by the square root of its weight,
$
\mathbf{P}=\begin{bmatrix}
\sqrt{w(p_1)}p_1 &
\cdots
\sqrt{w(p_n)}p_n
\end{bmatrix}^T
$
and let $\mathbf{p} = \sqrt{w(p)}p^T$ for every $p\in P$.
 The projection of the rows of $\mathbf{P}$ onto $S'$ is $\mathbf{P}X\in\ensuremath{\mathbb{R}}^{n\times k}$ using the column base of $X$, and $\mathbf{P}XX^T\in\ensuremath{\mathbb{R}}^{n\times d}$ in $\ensuremath{\mathbb{R}}^d$. Hence, for every $p\in P$, the weighted squared distance from $p$ to $S'$ is
\begin{align}
w(p)\cdot D(p,S') =  w(p) \norm{p^T-p^TXX^T}_2^2 = \norm{\sqrt{w(p)} \cdot p^T- \sqrt{w(p)} \cdot p^TXX^T}_2^2 =\norm{\mathbf{p} - \mathbf{p}XX^T}_2^2. \nonumber
\end{align}

By letting $Y\in \ensuremath{\mathbb{R}}^{d\times(d-k)}$ be the matrix that spans the orthogonal complement subspace of $S'$ (i.e., $Y^TY=I$ and $[X,Y][X,Y]^T=I$), we obtain
\begin{align*}
\mathbf{p} = \mathbf{p}I =\mathbf{p}(XX^T + YY^T) = \mathbf{p}XX^T + \mathbf{p}YY^T,
\end{align*}
and by subtracting $\mathbf{p}XX^T$ from both sides and applying squared norm we get that
\begin{align*}
w(p)D(p,S')=\norm{\mathbf{p}- \mathbf{p}XX^T}_2^2 = \norm{\mathbf{p}YY^T}^2_2 = \norm{\mathbf{p}Y}^2_2
\end{align*}
From the last equality we get that the sum of squared distance from the set $P$ to the $k$-subspace $S'$ is
\begin{align*}
\displaystyle{\sum_{q\in P}}  w(q)\cdot D(q,S') = \displaystyle{\sum_{q\in P}} \norm{\mathbf{q}Y}^2_2 =  \norm{\mathbf{P}Y}^2_F
\end{align*}

Thus by letting $\mathcal{S}_{d}$ be the set of all (non-affine) $k$-subspaces of $\ensuremath{\mathbb{R}}^d$, we get that the sensitivity of a point $p\in P$ in the query space $(P,w,\mathcal{S}_{d},D)$ is:
\begin{align}
s(p) = \sup_{S'\in \mathcal{S}_{d}} \frac{w(p)D(p,S')}{\displaystyle{\sum_{q\in P}} w(q)D(q,S')} =
 \sup_{Y\in \ensuremath{\mathbb{R}}^{d\times (d-k)},  Y^TY=I} \frac {\norm{\mathbf{p}Y}_2^2}{\norm{\mathbf{P}Y}_F^2},\label{sensreduction}
\end{align}
such that the denominator is not zero.

\newcommand{\textsc{Non-Affine-Sensitivity}}{\textsc{Non-Affine-Sensitivity}}
\setcounter{AlgoLine}{0}
\begin{algorithm}[th]
\caption{$\textsc{Non-Affine-Sensitivity}(P,w,p, k,\eps)$; see Lemma~\ref{lemma:svd}}\label{alg:tight1}
{\begin{minipage}{\textwidth}
\begin{tabbing}
\textbf{Input:} \quad\=A weighted set $(P,w)$ of $n$ points in $\mathbb{R}^d$, a point $p\in P$,
\\\> an integer $k\in\br{0,\cdots,d-1}$, and an error parameter $\eps\in(0,1)$. \\
\textbf{Output:} An additive $\eps$-approximation $s'$ to the sensitivity $s(p)$ of $p$
\end{tabbing}\end{minipage}}

\nl $\mathbf{p} :=  \sqrt{w(p)}{p}^{T}$  ; $\mathbf{P}:=\begin{bmatrix}
        \sqrt{w(p_1)}{p}_1  \cdots  \sqrt{w(p_n)}{p}_n \\
            \end{bmatrix}^T$ \label{matP}

\nl\If{$k=d-1$} {

\nl denote by $UDV^T$  the thin SVD of $\mathbf{P}.$

\nl\Return $\norm{u}_2^2$ \quad\quad\tcp{where $u$ is the corresponding row in $U$ for $\mathbf{p}$ in $\mathbf{P}$.} \label{u_i line bound}

}

\nl$\ell:= d-k$ \label{big p}

\nl$\displaystyle{\gamma := \frac{\sum_{i=1}^{\ell} \lambda_{d-i+1}( \mathbf{P}^T \mathbf{P}) }{\sum_{i=1}^{\ell} \lambda_{i}( \mathbf{P}^T \mathbf{P})}}$ \\ \tcp{ where $\lambda_{i}( \mathbf{P}^T \mathbf{P})$ is the $i$th eigenvalue of $ \mathbf{P}^T \mathbf{P}$ for every $i\in \br{1,\cdots,d}$.} \label{def gamma

$X :=$ any $d\times {\ell}$ matrix whose columns are orthonormal (i.e., $X^TX=I$).

$s_{new} :=\infty$ ; $s_{old} := -\infty $

\While {$\displaystyle {s_{new} -  \frac {\eps \gamma}{(1-\gamma)}\geq s_{old}}$}
{\label{while}

$\displaystyle {s_{new}:= \frac{\norm{{\mathbf{p}}X}^2_2}{\norm{{\mathbf{P}}X}^2_F}}$\label{snew}

$G := {\mathbf{p}^T \mathbf{p} } - s_{new} \cdot {\mathbf{P}^T \mathbf{P}}$\label{snewafter}

$X := \begin{bmatrix}
        x_1 \cdots  x_\ell \\
            \end{bmatrix}$  \quad\quad\tcp{a matrix that its cols are the eigenvectors corresponding to the $\ell$'th largest eigenvalues of $G$.}

}

$s' = s_{new}+\eps$

\Return $s'$

\end{algorithm}
\begin{lemma}\label{lemma:svd}
Let $(P,w)$ be a weighted set of $n$ points in $\ensuremath{\mathbb{R}}^d$, $k\in \br{0,\cdots,d-1}$ be an integer, $\eps \in (0,1)$ be an error parameter, and let $p\in P$. Let $\mathcal{S}_{d}$ denote the set of all non-affine $k$-subspaces in $\ensuremath{\mathbb{R}}^d$,  $\displaystyle{s(p) = \sup_{S'\in \mathcal{S}_{d}} \frac{w(p)D(p,S')}{\displaystyle{\sum_{q\in P}} w(q)D(q,S')}}$ denote the sensitivity of $p$ in the $(P,w,\mathcal{S}_{d},D)$ query space, where the denominator is not zero, and let $s'$ be the output of a call to $\textsc{Non-Affine-Sensitivity}(P,w,p, k,\eps)$; See Algorithm~\ref{alg:tight1}. Then the following holds according to $k$:
\begin{enumerate}[(i)]
\item if $k=d-1$ then  $s' = s(p)$.
\item if $k\in \br{0,\cdots,d-2}$ then  $s(p) \leq s' \leq s(p) + \eps$.
\end{enumerate}
\end{lemma}

\begin{proof}
Let $\mathbf{P}$ and $\mathbf{p}$ be the matrix and row vector that are defined at Line~\ref{matP} of Algorithm~\ref{alg:tight1}.

\textbf{Proof of Claim (i) $k=d-1$}:
 Let $\mathbf{P}=UDV^T$ denote the thin Singular Value Decomposition of $\mathbf{P}$,
 and let $\mathbf{p}=uD V^T$ (where $u$ is the corresponding row in $U$ for the row $\mathbf{p}$ in $\mathbf{P}$ ). Let $x^*$ be the vector that maximizes $\frac{\abs{\mathbf{p}x}^2}{\norm{\mathbf{P}x}_2^2}$ over every $x\in\ensuremath{\mathbb{R}}^d$ such that $\norm{\mathbf{P}x}_2>0$ and $x^Tx=1$. It is well known (e.g.~\cite{yang2017weighted}) that,\\
\begin{align*}
s(p) =
&\frac{\abs{\mathbf{p}x^*}^2}{\norm{\mathbf{P}x^*}_2^2}
= \frac{\abs{uD V^Tx^*}^2}{\norm{UD V^Tx^*}_2^2}
= \frac{\abs{uDV^Tx^*}^2}{\norm{DV^Tx^*}_2^2}\leq \frac{\norm{u}_2^2\norm{DV^Tx^*}_2^2}{\norm{DV^Tx^*}_2^2}
= \norm{u}_2^2,
\end{align*}
where the first equality holds by~\eqref{sensreduction}, the second is by the definition of $u,U,D$ and $V$, the third is since the columns of $U$ are orthonormal, and the inequality holds by the Cauchy Schwarz inequality.
We now prove that this upper bound is tight. Indeed, substituting $x=VD^{-1}u^T$ (where $D^{-1}$ is the inverse matrix of $D$) attains this maximum as,
\begin{align*}
\frac{\abs{\mathbf{p}x}^2}{\norm{\mathbf{P}x}_2^2} =
\frac{\abs{uDV^TVD^{-1}u^T}^2}{\norm {DV^TVD^{-1}u^T}_2^2}=
\frac{\abs{uu^T}^2}{\norm{u^T}_2^2} = \norm{u}_2^2.
\end{align*}
Thus, $s(p) =\norm{u}_2^2 $, which is the returned value of Algorithm~\ref{alg:tight1} for the case $k=d-1$.; See Line~\ref{u_i line bound}.
\newcommand{\X}{\displaystyle{ X^*}}\\

\textbf{Proof of Claim (ii) $k\in\br{0,\cdots,d-2}$}: Let $\X$ be the Matrix that maximizes $\frac{\norm{\mathbf{p}X}_2^2}{\norm{\mathbf{P}X}_F^2}$ over every $X\in\ensuremath{\mathbb{R}}^{d\times (d-k)}$ such that $\norm{\mathbf{P}X}_F>0$ and $X^TX=I$. We have,
\begin{align*
s(p) =
\frac{\norm{\mathbf{p}\X}_2^2 }{ \norm{\mathbf{P}{\X}}_F^2}=
{\frac{Tr({\X}^T\mathbf{p}^T\mathbf{p}{\X})}{Tr({\X}^T \mathbf{P}^T \mathbf{P} {\X})}}
\end{align*}

Let $j_{\max}>1$ be the number of iterations that are executed in the ``while" loop of Algorithm~\ref{alg:tight1} until it stops, and let $s_j$ be the value of $s_{new}$ during the execution of Line~\ref{snewafter} in the $j$th iteration for every $j\in \br{1,\cdots,j_{\max}}$.
The while loop of Algorithm~\ref{alg:tight1} is the same while loop of Algorithm~2 in~\cite{zhang2010fast}, where the main difference is the stopping criterion.
  In \cite{zhang2010fast} it was proven that $s_{new}$ (in the ``while`` loop of Algorithm~\ref{alg:tight1}) converges to the global supremum of $$ \displaystyle{\frac{Tr(X^T\mathbf{p}^T\mathbf{p}X)}{Tr(X^T {\mathbf{P}^T\mathbf{P}} X)}}.$$
 Moreover $s_j \leq s(p)$ for every $j\in \br{1,\cdots,j_{\max}}$.

  Let $\gamma$ be defined as at Line~\ref{def gamma} of Algorithm~\ref{alg:tight1}. Based on Theorem 5.1 in \cite{zhang2010fast}, for any integer $j\geq 1$, we have that,
\begin{align*}
s(p) - s_j &\leq (1-\gamma)(s(p)  - s_{j-1}) .
\end{align*}
Rearranging yields,
\begin{align}
s(p) \leq \frac{s_j - s_{j-1}(1-\gamma)}{\gamma}. \label{th5.1ofzhang}
\end{align}
By the stopping criterion (in Line~\ref{while} of Algorithm~\ref{alg:tight1}) we have that after the last iteration
\begin{align}\left(s_{j_{\max}} -  \frac { \varepsilon \gamma}{1-\gamma}\right)< s_{{j_{\max}-1}}. \label{sc-1 bound}\end{align}
By combining~\eqref{sc-1 bound} with~\eqref{th5.1ofzhang} we obtain an upper bound on $s(p)$, as:
$$ s(p) \leq  \frac{s_{j_{\max}}  -  (s_{j_{\max}}  -  \frac {  \eps \gamma}{1-\gamma})  (1-\gamma)}{\gamma} = \frac{s_{j_{\max}}  - s_{j_{\max}}  + \gamma s_{j_{\max}}  +  \eps\gamma}{\gamma}=s_{j_{\max}} + \eps. $$
By the above inequality and since $s_{j_{\max}} \leq s(p)$, we have
\begin{align*}
s(p) \leq  s_{j_{\max}}  +  \eps \leq s(p)+ \eps.
\end{align*}
We conclude that the returned value of the algorithm $s' =s_{j_{\max}} +  \eps$ satisfies Claim~(ii) as,
$$s(p) \leq  s' \leq s(p)+ \eps.$$
\end{proof}
\section{Reduction from Affine to Non-Affine Subspace} \label{reducsec}
In this section we use the two integers $n,d\geq1$,  the number $z\geq1$ and $\eps\in(0,  \frac{1}{2^{z+1} +2} ]$ as the additive error of our sensitivity bounds, which is polynomial in $1/n$ in our experimental results. We let $(P,w)$ be a weighted set of $n$ points in $\ensuremath{\mathbb{R}}^d$, $\psi= \left(\frac{\eps}{z}\right)^z$,  $\displaystyle r = 1+ \max_{p\in P}\frac{D_z(p,0_d)}{\psi\eps^2}$, and $\displaystyle{e_{d+1}=(0, \cdots,0,1)\in\ensuremath{\mathbb{R}}^{d+1}}$.

For every $q\in P$ we let $q' = (q\mid r) $ and $P' = \br{ q' \mid q\in P }$. The set of all affine $k$-subspaces in $\ensuremath{\mathbb{R}}^d$ is denoted by $\mathcal{S}^{A}_{d}$.
  For every affine $k$-subspace $S\in \mathcal{S}^{A}_{d}$, we define $\displaystyle{S'' := \{ (x\mid r)\mid x\in S\}}$ to be its corresponding affine $k$-subspace in $\ensuremath{\mathbb{R}}^{d+1}$, and $S'$ to be the corresponding (non-affine) $(k + 1)$-subspace of $\ensuremath{\mathbb{R}}^{d+1}$ that is spanned by $S''$; see Figs.~\ref{fig:big} and~\ref{fig:small}. Finally let $\displaystyle{\mathcal{S}_{d+1}}$ denote the union over all $(k+1)$ non-affine subspaces of $\ensuremath{\mathbb{R}}^{d+1}$.

\begin{figure}
  \centering
  \includegraphics[width=0.8\columnwidth]{big.png}
  \caption{ Illustration of Claim (i) at Lemma~\ref{theorem:j-affine}.  This figure illustrates the case where $d=1$, $z=1$, and $k=0$, $p$ is a point in $\ensuremath{\mathbb{R}}$ such that $D_1(p,0_d) = \mathrm{dist}(p,0_d) = \abs{p-0_{d}} \leq \eps^2 r$, and $p'=(p\mid r)\in \ensuremath{\mathbb{R}}^2$, $S$ is an affine $0$-subspace of $\ensuremath{\mathbb{R}}$ (a point in $\ensuremath{\mathbb{R}}$ for this case) such that $D_1(0_d,S)  = \mathrm{dist}(0_d,S) = \abs{S-0_{d}}\geq \eps r$ , $S''$ is the corresponding affine $0$-subspace of $\ensuremath{\mathbb{R}}^2$, and $S'$ is the (non-affine) $1$-subspace (i.e., line)  of $\ensuremath{\mathbb{R}}^2$ that passes through $S''$. This illustration aims to show that the distance from $p$ to $S$ (i.e, $D_1(p,S)$) is approximately equal to the distance from $0_d$ to $S$ ($D_1(p,S)$), and that $D_1(p',S')$ is approximately equal to $D_1(re_{2},S' )$.  }
  \label{fig:big}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.8\columnwidth]{small.png}
  \caption{Illustration of Claim (ii) at Lemma~\ref{theorem:j-affine}. Same as Figure~\ref{fig:big} we have $d=1$, $z=1$, $k=0$, $p$ is a point in $\ensuremath{\mathbb{R}}$ such that $D_1(p,0_d) = \mathrm{dist}(p,0_d) = \abs{p-0_{d}} \leq \eps^2 r$, $p'=(p\mid r)\in \ensuremath{\mathbb{R}}^2$, $S$ is an affine $0$-subspace of $\ensuremath{\mathbb{R}}$ (a point in $\ensuremath{\mathbb{R}}$ for this case) such that $D_1(0_d,S)  = \mathrm{dist}(0_d,S) = \abs{S-0_{d}}< \eps r$ , $S''$ is the corresponding affine $0$-subspace of $\ensuremath{\mathbb{R}}^2$, and $S'$ is the (non-affie) $1$-subspace (i.e., line)  of $\ensuremath{\mathbb{R}}^2$ that passes through $S''$. This illustration aims to show that the distance from $p$ to $S$ (i.e, $D_1(p,S)$) is approximately equal to the distance from $p'$ to $S'$ (i.e., $D_1(p',S')$).}
\label{fig:small}
\end{figure}
\begin{lemma}\label{theorem:j-affine}
Let $S\in \mathcal{S}^{A}_{d}$ be an affine $k$-subspace of $\ensuremath{\mathbb{R}}^d$. For every $p \in P$ and its corresponding point $p'=(p\mid r)$, Claims (i)--(ii) hold as follows:
\begin{enumerate}[(i)]
\item if $D_z(0_d,S) \geq \eps r$ then
\begin{align}
&| D_z(0_{d},S) - D_z(p,S) | \leq 2\eps \cdot D_z(0_{d},S)\label{oo1}
\intertext{and}
&| D_z(re_{d+1},S') - D_z(p',S') |\leq  (2^z+1 )\eps \cdot D_z(re_{d+1},S').\label{oo2}
\end{align}
\item if $D_z(0_d,S) < \eps r$, then
\begin{align}
&D_z(p',S') \leq D_z(p,S) \leq (1+\constalaa \eps)D_z(p',S').\label{oo3}
\end{align}
\end{enumerate}
\end{lemma}

\begin{proof}
By~\cite[Lemma~2.1]{Feldman:2012:DRW:2095116.2095222}, for every $x,y\in\ensuremath{\mathbb{R}}^d$ we have
\begin{align}
\abs{D_z(x,S) - D_z(y,S)} \leq \frac{D_z(x,y)}{\psi} + \eps D_z(x,S).\label{knownfact}
\end{align}
Let $p\in P$ and let $p'=(p\mid r)$.

\textbf{Proof of Claim (i): } Inequality~\eqref{oo1} holds since
\begin{align*}
|D_z(p,S) - D_z(0_d,S)|&\leq \frac{D_z(p,0_d)}{\psi} + \eps  D_z(0_d,S) \leq \eps^2r + \eps  D_z(0_d,S) \leq 2\eps D_z(0_d,S),
\end{align*}
where the first inequality holds by substituting $x=0_d$ and $y=p$ in~\eqref{knownfact}, the second is by the definition of $r$, and the last inequality holds by the assumption of Claim~(i).

Let $h_{S''}=\mathrm{proj}(re_{d+1},S'')$. To prove~\eqref{oo2}, we observe that
\begin{align}
\mathrm{dist}(re_{d+1}, h_{S''}) = \mathrm{dist}(re_{d+1}, S'') = \mathrm{dist}(0_d,S)=D^{1/z}_z(0_d,S) \geq (\eps r)^{\frac{1}{z}}, \label{dist bound}
\end{align}
where the first equality holds by the definition of $h_{S''}$, the second holds by the definition of $S''$ and $e_{d+1}$, the third holds by the definition of $D_z$, and the inequality holds by taking the power of $1/z$ of each side of the assumption $D_z(0_d,S)\geq \varepsilon r$.
In addition let $h_{S'}=\mathrm{proj}(re_{d+1},S')$, and let $\beta_{S'} \in [0, \pi/2)$ denote the angle $\angle(re_{d+1}, 0_{d+1}, h_{S'})$. Hence
\begin{align}
& \mathrm{dist}(re_{d+1},h_{S'}) = \sin{\beta_{S'}} \cdot \mathrm{dist}(re_{d+1}, 0_{d+1}) = \sin{\beta_{S'}} \cdot r.\label{dist(re1,projjre)}
\end{align}
Now we compute a lower bound on $\sin{\beta_{S'}}$:
\begin{align}
\sin{\beta_{S'}} & =  \frac{\mathrm{dist}(re_{d+1},h_{S''} )} {\mathrm{dist}(0_{d+1} , h_{S''})} \label{eq1}
\\&\geq  \frac{\mathrm{dist}(re_{d+1}, h_{S''})} { \mathrm{dist}( 0_{d+1},re_{d+1})+  \mathrm{dist}(re_{d+1}, h_{S''})  }\label{eq2}
\\&\geq  \frac{(\eps r)^{\frac{1}{z}}} { \mathrm{dist}( 0_{d+1},re_{d+1})+  (\eps r)^{\frac{1}{z}} }\label{eq3}
\\&=  \frac{(\eps r)^{\frac{1}{z}}}{r+ (\eps r)^{\frac{1}{z}}}
=\frac{\eps^{\frac{1}{z}}} {\frac{r}{r^{\frac{1}{z}}} +\eps^{\frac{1}{z}}}
=\frac{\eps^{\frac{1}{z}}} {r^{\frac{z-1}{z}} +\eps^{\frac{1}{z}}}\label{minval}\\
&\geq \frac{\eps^{\frac{1}{z}}}{2r^\frac{z-1}{z}} \label{minval1}
\end{align}
where~\eqref{eq2} holds by the triangle inequality, \eqref{eq3} holds by~\eqref{dist bound},

and~\eqref{minval} holds since $\mathrm{dist}(0_{d+1},re_{d+1})= r$.
 Equation \eqref{minval1} holds since $r^{\frac{z-1}{z}}\geq\eps^{\frac{1}{z}}$ for $z\geq1$, $r>1$ and $\eps<1$.
 Plugging~\eqref{minval1} in~\eqref{dist(re1,projjre)} yields,
\begin{align}
& \mathrm{dist}(re_{d+1},h_{S'}) \geq  \frac{\eps^{\frac{1}{z}}}{2r^\frac{z-1}{z}}\cdot r = \frac{(\eps r)^{1/z}}{2}.\label{distbound}
\end{align}
By~\eqref{distbound} and the definition of $D_z$,
\begin{align}
& D_z(re_{d+1},S') = D_z(re_{d+1},\mathrm{proj}(re_{d+1},S'))=D_z(re_{d+1},h_{S'})\geq\bigg(\frac{(\eps r)^{1/z}}{2}\bigg)^z
=\frac{\eps r}{2^z}.\label{goodone}
\end{align}
To complete the proof of Claim (i),
\begin{align}
|D_z(p',S') - D_z(re_{d+1},S')| &\leq \frac{D_z(p',re_{d+1})}{\psi} + \eps D_z(re_{d+1},S') \label{nn1} \\
& = \frac{D_z(p ,0_{d})}{\psi} + \eps D_z(re_{d+1},S')\label{nn2} \\
& \leq\eps^2 r + \eps  D_z(re_{d+1},S') \label{nn3} \\
& \leq   (2^z+1 )\eps D_z(re_{d+1},S'), \label{diffbound2}
\end{align}
where~\eqref{nn1} holds by substituting $x=p'$, $y=re_{d+1}$ and $S=S'$ in~\eqref{knownfact},~\eqref{nn2} is by the definition of $e_{d+1}$ and $p'$,~\eqref{nn3} is by the definition of $r$, and~\eqref{diffbound2} holds by~\eqref{goodone}.
\newcommand{ h_{S''}^{p'}}{ h_{S''}^{p'}}
\newcommand{\alpha_{S''}^{p'}}{\alpha_{S''}^{p'}}

\textbf{Proof of Claim (ii):} In this case $D_z(0_d,S)< \eps r$. Let $h_{S''}^{p'}=\mathrm{proj}(p',S'')=(\mathrm{proj}(p,S)^T\mid r)^T$ and $\alpha_{S''}^{p'} \in [0, \pi/2)$ denote the angle $\angle(\mathrm{proj}(p', S'), p', h_{S''}^{p'})$. Hence,
\begin{align}\label{p,S}
\mathrm{dist}(p,S) = \mathrm{dist}(p',S'')=\frac{\mathrm{dist}(p',S') }{\cos{\alpha_{S''}^{p'}}}  .
\end{align}
We have that $\cos{\alpha_{S''}^{p'}} \leq 1$, thus $\mathrm{dist}(p',S') \leq \mathrm{dist}(p,S)$. Observing that $$\alpha_{S''}^{p'}=\angle(re_{d+1}, 0_{d+1}, h_{S''}^{p'}) =\angle(re_{d+1}, 0_{d+1}, h_{S'}) = \beta_{S'}$$ yields,
\begin{align}
\cos{\alpha_{S''}^{p'}} =  \frac{re_{d+1}^T  h_{S''}^{p'}}{\norm{re_{d+1}}_2\norm{ h_{S''}^{p'}}_2} = \frac{(0,\cdots,0,r) (\mathrm{proj}(p,S)^T\mid r)^T}{\norm{(0,\cdots,0,r)}_2\norm{ h_{S''}^{p'}}_2} = \frac{r^2}{r\norm{ h_{S''}^{p'}}_2} = \frac{r}{\norm{ h_{S''}^{p'}}_2}.\label{cosal}
\end{align}
By the triangle inequality,
\begin{align}\label{normh}
\norm{ h_{S''}^{p'}}_2 \leq \norm{0_{d+1} -re_{d+1}}_2 + \norm{re_{d+1} - p'}_2 + \norm{p' -  h_{S''}^{p'}}_2.
\end{align}
We also have $\norm{p'- h_{S''}^{p'}}_2=\mathrm{dist}(p,S)$  and by the triangle inequality $\mathrm{dist}(p,S)\leq \norm{p}_2 + \mathrm{dist}(0_d,S)$. The bound on $\norm{p}_2$ is,
\begin{align}
\norm{p}_2 = (D(p,0_d))^{1/z} \leq (\eps^2\psi r)^{1/z}\leq \eps r. \label{normp}
\end{align}
Thus,
\begin{align}\label{boundph}
\norm{p'- h_{S''}^{p'}}_2=\mathrm{dist}(p,S)\leq \norm{p}_2 + \mathrm{dist}(0_d,S) \leq 2\eps r,
\end{align}
where the last inequality is by~\eqref{normp} and the assumption of Claim (ii). And by plugging~\eqref{boundph} in~\eqref{normh} and using the facts that $\norm{0_{d+1} -re_{d+1}}_2=r$ and $\norm{re_{d+1} - p'}_2 = \norm{p}_2\leq \eps r$ we get that,
\begin{align*}
\norm{ h_{S''}^{p'}}_2 \leq r +\norm{p}_2 + 2\eps r \leq (1+3\eps) r.
\end{align*}
Together with~\eqref{cosal} this yields a lower bound on $\cos{\alpha_{S''}^{p'}}$ as,
\begin{align*}
\cos{\alpha_{S''}^{p'}}= \frac{r}{\norm{ h_{S''}^{p'}}_2} \geq  \frac{1}{ 1+3\eps }.
\end{align*}
Now we obtain an upper bound for~\eqref{p,S} as,
\begin{align*}
\mathrm{dist}(p,S) = \frac{\mathrm{dist}(p',S') }{\cos{\alpha_{S''}^{p'}}}  \leq (1+3 \eps) \mathrm{dist}(p',S') .
\end{align*}

Thus we have that $ \mathrm{dist}(p',S') \leq \mathrm{dist}(p,S)  \leq(1+3 \eps) \mathrm{dist}(p',S')$. By the definition of $D_z$ we have
\begin{align}
D_z(p',S') \leq D_z(p,S) \leq (1+3\eps)^z D_z(p',S'). \label{plughere}
\end{align}
This proves Claim (ii) for the case $z=1$. Otherwise, by the Bernoulli's inequality we have that $(1+x)^y \leq 1+ \frac{xy}{1-(y-1)x}$ for every $x\in [-1,\frac{1}{y-1})$ and $y>1$. Thus substituting $y=z,x=3\eps$ yields,
\begin{align*}
(1+3\eps)^z D_z(p',S')  \leq \left(1 + \frac{3z\eps}{1-(z-1)3\eps}\right) D_z(p',S').
\end{align*}
Observe that since $\eps \leq  \frac{1}{2^{z+1} +2} \leq \frac{1}{6(z-1)}$ we have that $1-(z-1)3\eps \geq 1-(z-1)\frac{3}{6(z-1)} = 1/2$.  Hence,
\begin{align*}
 (1+3\eps)^z D_z(p',S')  \leq \left(1 + \frac{3z\eps}{1-(z-1)3\eps}\right) D_z(p',S') \leq (1 +\constalaa \eps) D_z(p',S').
\end{align*}
Plugging the last inequality in~\eqref{plughere} proves Claim (ii) for $z>1$ as
$$D_z(p',S') \leq D_z(p,S) \leq  (1 +\constalaa \eps) D_z(p',S').$$
\end{proof}
\begin{lemma}\label{affinetononaffinecor}
For every $\displaystyle{p\in P}$ and its corresponding point $p'=(p\mid r)$ we hav
\begin{align*}
\left|\sup_{S\in \mathcal{S}^{A}_{d}}\frac{w(p)\cdot D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)} - \sup_{S'\in \mathcal{S}_{d+1}}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}\right| \leq 16\eps(2^z +1) \cdot \sup_{S'\in \mathcal{S}_{d+1}}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}.
\end{align*}
\end{lemma}

\begin{proof}
For every (non-affine) subspace $S'\in \mathcal{S}_{d+1}$ we denote $h_{S'}=\mathrm{proj}(re_{d+1},S')$ and $\beta_{S'} \in [0, \pi/2]$ denote the angle $\angle(re_{d+1}, 0_{d+1}, h_{S'})$. The following observation is from Figure~\ref{fig:big}:\\
\begin{observation} \label{obs:1} Let $S'$ be a non-affine $(k+1)$-subspace of $\ensuremath{\mathbb{R}}^{d+1}$ such that $\beta_{S'} < \pi/2$. Then the intersection of $S'$ with the hyperplane $\br{(x\mid r) \mid x \in \ensuremath{\mathbb{R}}^d}$ is the affine $k$-subspace $S'' \subset S'$ of $\ensuremath{\mathbb{R}}^{d+1}$ such that $S'$ is the linear span of $S''$, and the $(d+1)$th (last) coordinate of every point $x\in S''$ is $r$. Moreover there is an affine $k$-subspace $S$ of $\ensuremath{\mathbb{R}}^d$ such that $S'' = \br{(x \mid r) \mid x \in S}$. Hence,
\begin{align}
\sin{\beta_{S'}} = \frac{\mathrm{dist}(re_{d+1} , S'')}{ \mathrm{dist}(0_{d+1} , S'') }= \frac{\mathrm{dist}(re_{d+1} , S'')}{\sqrt{\mathrm{dist}^2(re_{d+1} , S'') + \mathrm{dist}^2(re_{d+1},0_{d+1})^2 } } = \frac{\mathrm{dist}(0_{d} , S)}{\sqrt{ \mathrm{dist}^2(0_{d} , S) +r^2}}. \label{sineq}
\end{align}
\end{observation}
The above observation will be used through the proof. Let $c_0 =\frac{(\eps r)^{1/z}}{ \sqrt{(\eps r)^{2/z} + r^2}}$. We partition the query set $\mathcal{S}_{d+1}$ into two disjoint subsets:
\begin{enumerate}[(i)]
\item $\displaystyle {Q'_{0} = \br{S' \in \mathcal{S}_{d+1} \mid  \sin{\beta_{S'}} \geq  c_0  }}$, and
\item $\displaystyle {Q'_{1} = \br{S' \in \mathcal{S}_{d+1}\mid 0\leq\sin{\beta_{S'}} < c_0  }}$.
\end{enumerate}
Similarly, we partition $\mathcal{S}^{A}_{d}$ into :
\begin{enumerate}[(i)]
\item $\displaystyle{Q_{0} =\br{S \in \mathcal{S}^{A}_{d}  \mid  D_z(0_d,S) \geq \eps r}}$, and
\item $\displaystyle{Q_{1} =\br{S \in \mathcal{S}^{A}_{d}  \mid  D_z(0_d,S) < \eps r }}$.
\end{enumerate}
Let $c_1 = 2^z +1$. We first proof the following pair of claims: \\
\textbf{Claim (i).} For every $p\in P$ we have
\begin{align*}
\left|\sup_{S\in Q_0}\frac{w(p)\cdot D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)} - \sup_{S'\in Q_0'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}\right| \leq 16c_1\eps \sup_{S'\in Q_0'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}.
\end{align*}
\textbf{Claim (ii).} For every $p\in P$ we have
\begin{align*}
\left|\sup_{S\in Q_1}\frac{w(p)\cdot D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)} - \sup_{S'\in Q_1'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}\right| \leq 16c_1\eps \sup_{S'\in Q_1'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}.
\end{align*}
Since $\mathcal{S}^{A}_{d} = Q_0 \cup Q_1$ and $\mathcal{S}_{d+1} = Q_0' \cup Q_1'$, combining both claims poofs the lemma.

\textbf{Proof of Claim (i).} By~\eqref{oo1} and since $c_1\geq 2$, for every pair $S\in Q_0$ and $p\in P$ we have
$$\displaystyle {(1-c_1\eps)D_z(0_{d},S) \leq D_z(p,S) \leq(1+c_1\eps) D_z(0_{d},S)}.$$
Hence,
\begin{align*}
\frac{w(p) (1-c_1\eps)D_z(0_d,S) }{\displaystyle{\sum_{q\in P}} w(q)  (1+c_1\eps) D_z(0_d,S) } \leq \frac{w(p)D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)}  \leq \frac{w(p)(1+ c_1\eps)D(0_d,S) }{\displaystyle{\sum_{q\in P}} w(q) (1- c_1\eps)D_z(0_d,S) }.
\end{align*}
Since the above inequality holds for every $S \in Q_0$ we get that,
\begin{align}
\frac{(1- c_1\eps)}{(1+c_1\eps)}\frac{w(p)}{\displaystyle{\sum_{q\in P}} w(q) } \leq \sup_{S\in Q_0}\frac{w(p)D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)}  \leq \frac{(1+c_1\eps)}{(1- c_1\eps)}\frac{w(p)}{ \displaystyle{\sum_{q\in P}} w(q) }. \label{beforebound}
\end{align}

Let $p'\in P'$. We now prove that for every $S' \in Q_0'$
\begin{align}
(1-c_1\eps)D_z(re_{d+1},S') \leq D_z(p',S') \leq(1+c_1\eps) D_z(re_{d+1},S') \label{bigsin}
\end{align}
by case analysis: first for $\sin\beta_{S'} \in [c_0,1)$ and then for $\sin\beta_{S'}=1$.
 If $S' \in  \br{ S' \in \mathcal{S}_{d+1} \mid  \sin{\beta_{S'}}\in [c_0,1) }$ then we have that $\beta_{S'} < \pi/2$. Hence, by  Observation~\ref{obs:1} there is an affine $k$-subspace $S$ (of $\ensuremath{\mathbb{R}}^d$) such that
$$\displaystyle{\sin{\beta_{S'}} = \frac{\mathrm{dist}(0_{d} , S)}{\sqrt{ \mathrm{dist}^2(0_{d} , S) + r^2}}}.$$
Combining this equality with the fact that $\sin{\beta_{S'}} \geq  c_0 =\frac{(\eps r)^{1/z}}{\sqrt{(\eps r)^{2/z} + r^2}}$ yields that $\mathrm{dist}(0_{d} , S) \geq (\eps r)^{1/z}$. Taking the power of $z$ from both sides yields $D_z(0_d,S) \geq \eps r$. Using this in~\eqref{oo2} yields that~\eqref{bigsin} holds for the case $\sin\beta_{S'} \in [c_0,1)$.

If $\sin \beta_{S'}= 1$ then we have that $\beta_{S'} = \pi/2$. This implies that $h_{S'} = 0_{d+1}$. Hence,
\begin{align}
\mathrm{dist}(re_{d+1},S') = \mathrm{dist}(re_{d+1},h_{S'}) =\mathrm{dist}(re_{d+1},0_{d+1}) = r \label{eqr}.
\end{align}
Hence, for every $S' \in \br{ S' \in\mathcal{S}_{d+1} \mid \sin{\beta_{S'}}  = 1 }$ we have
\begin{align}
\abs{D_z(re_{d+1},S') - D_z(p',S')}
&\leq \frac{D_z(re_{d+1},p')}{\psi} + \eps D_z(re_{d+1},S')\label{n0}   \\
&\leq \frac{D_z(0_{d},p)}{\psi} + \eps D_z(re_{d+1},S')\label{n1} \\
&\leq \frac{r\psi\eps^2}{\psi} +  \eps D_z(re_{d+1},S')\label{n2}\\
&= \eps^2 r + \eps  D_z(re_{d+1},S') \nonumber \\
&\leq  2\eps  D_z(re_{d+1},S') \label{n4} \\
& \leq  c_1\eps  D_z(re_{d+1},S')\label{n5} ,
\end{align}
where~\eqref{n0} holds by substituting $x=re_{d+1}$, $y=p'$ and $S=S'$ in~\eqref{knownfact},~\eqref{n1} is by the definition of $r$ and $p'$,~\eqref{n2} holds by the definiton of $\psi$ and $r$, ~\eqref{n4} holds by~\eqref{eqr}, and ~\eqref{n5} holds since $c_1 > 2$. This proves~\eqref{bigsin} also for the case that $\sin \beta_{S'}=1$.  Hence,~\eqref{bigsin} holds for every $S' \in Q_0'$.

By~\eqref{bigsin} we get that for every $p'\in P'$
\begin{align}
\frac{(1-c_1\eps)}{(1+ c_1\eps)}\frac{w(p)  }{\displaystyle{\sum_{q\in P}} w(q) } \leq \sup_{S'\in Q_0'}\frac{w(p) D_z(p',S') }{\displaystyle{\sum_{q\in P}} w(q)D_z(q',S') }  \leq \frac{(1+c_1\eps)}{(1- c_1\eps)}\frac{w(p)}{\displaystyle{\sum_{q\in P}} w(q)  } .\label{doneHtag}
\end{align}

Integrating~\eqref{doneHtag} with~\eqref{beforebound} yields
\begin{align}
\frac{(1- c_1\eps)^2}{(1+ c_1\eps)^2}\sup_{S'\in Q_0'} \frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')} &  \leq
\frac{(1- c_1\eps)}{(1+ c_1\eps)}\frac{w(p)}{\displaystyle{\sum_{q\in P}} w(q) } \label{re1}
\\ &  \leq \sup_{S\in Q_0}\frac{w(p)D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)} \label{re2} \\
&\leq \frac{(1+c_1\eps)}{(1- c_1\eps)}\frac{w(p)}{ \displaystyle{\sum_{q\in P}} w(q)  }\label{re3} \\
&\leq \frac{(1+ c_1\eps)^2}{(1- c_1\eps)^2}\sup_{S'\in Q_0'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')},  \label{re4}
\end{align}
where~\eqref{re1} holds by multiplying the right hand side of~\eqref{doneHtag} by $\frac{(1- c_1\eps)^2}{(1+ c_1\eps)^2}$,~\eqref{re2} and~\eqref{re3} hold by~\eqref{beforebound}, and~\eqref{re4} holds by multiplying the left hand side of~\eqref{doneHtag} by $\frac{(1+ c_1\eps)^2}{(1- c_1\eps)^2}$.
We have
\begin{align}
\frac{(1-c_1\eps)^2}{(1+ c_1\eps)^2} &= \frac{1+c_1^2\eps^2-2c_1\eps}{(1+c_1\eps)^2} =\frac{1+2c_1\eps+c_1^2\eps^2-4c_1\eps}{(1+c_1\eps)^2}\label{vav}
\\&=1-\frac{4c_1\eps}{(1+c_1\eps)^2} \geq  1-4c_1\eps \geq 1-16c_1\eps, \nonumber
\end{align}
and
\begin{align}
\frac{(1+c_1\eps)^2}{(1 -  c_1\eps)^2}&=\frac{1+c_1^2\eps^2+2c_1\eps}{(1 -c_1\eps)^2}=\frac{1 -2c_1\eps+c_1^2\eps^2 + 4c_1\eps}{(1 -c_1\eps)^2}\label{vav1}  \\&=1+\frac{4c_1\eps}{(1 -c_1\eps)^2} \leq 1+\frac{4c_1\eps}{(1 -\frac {c_1}{2 c_1} )^2} = 1+\frac{4c_1\eps}{ \frac {1}{4} } = 1+16c_1 \eps, \nonumber
\end{align}
where the inequality in~\eqref{vav1} holds since $\eps \in \left(0,\frac{1}{2(2^z +1)}\right)$. By plugging~\eqref{vav} and~\eqref{vav1} in \eqref{re4} we get
\begin{align*}
(1- 16c_1\eps)\sup_{S'\in Q_0'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')} \leq  \sup_{S\in Q_0} \frac{w(p)D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)}  \leq
(1+16c_1\eps)\sup_{S'\in Q_0'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}.
\end{align*}
This proves Claim (i) as
$$\left|\sup_{S\in Q_0}\frac{w(p)\cdot D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)} - \sup_{S'\in Q_0'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}\right| \leq 16c_1\eps \sup_{S'\in Q_0'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}.$$

\textbf{Proof of Claim (ii).}
Let $S'\in Q'_1$. We have that $\beta_{S'}<\pi/2$. Hence, by Observation~\ref{obs:1} there is an affine $k$-subspace $S$ (of $\ensuremath{\mathbb{R}}^d$) such that
\begin{align}
\displaystyle{\sin{\beta_{S'}} = \frac{\mathrm{dist}(0_{d} , S)}{\sqrt{ \mathrm{dist}^2(0_{d} , S) + r^2}}}. \label{sinbetaa}
\end{align}
By the definition of $S'$ we have that $\sin{\beta_{S'}} <c_0 = \frac{(\eps r)^{1/z}}{ \sqrt{(\eps r)^{2/z} +  r^2}}$. Combining this with~\eqref{sinbetaa} yields that $\mathrm{dist}(0_{d} , S) < (\eps r)^{1/z}$. Taking the power of $z$ from both sides yields that $D_z(0_d,S) < \eps r$.
Using this is~\eqref{oo3} yields that for every $p' \in P$ and $S'\in Q_1'$ we have
\begin{align*}
\frac{1}{(1+\constalaa\eps)}\frac{ w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')} \leq  \frac{w(p) D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)} \leq  (1+\constalaa\eps)\frac{w(p)  D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S') }.
\end{align*}
By combining the last inequality with the fact that
\begin{align*}
1-\constalaa\eps\leq1-\frac{\constalaa\eps}{1+\constalaa\eps}=\frac{1+\constalaa\eps-\constalaa\eps}{1+\constalaa\eps}=\frac{1}{1+\constalaa\eps},
\end{align*}
we obtain that
\begin{align*}
&(1-\constalaa\eps)\frac{ w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')} \leq  \frac{w(p) D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)} \leq  (1+\constalaa\eps)\frac{w(p)  D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}.
\end{align*}
Since the above inequalities hold for every $S' \in Q_1'$ and $16c_1=16(2^z+1)>\constalaa$, this proves Claim~(ii) as
\begin{align*}
(1-16c_1 \eps)\sup_{S' \in Q_1'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')} \leq\sup_{S \in Q_1} \frac{w(p)D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)}  \leq
(1+16c_1\eps)\sup_{S' \in Q_1'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')},
\end{align*}
i.e.,
\begin{align*}
\bigg|\sup_{S'\in Q_1'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')} -  \sup_{S\in Q_1}\frac{w(p)D_z(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D_z(q,S)}\bigg|  \leq 16\eps c_1 \cdot \sup_{S'\in Q_1'}\frac{w(p)D_z(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D_z(q',S')}.
\end{align*}
\end{proof}

\section{Sensitivity of affine $k$-subspaces}
\newcommand{\textsc{Affine-Sensitivity}}{\textsc{Affine-Sensitivity}}
\setcounter{AlgoLine}{0}
\begin{algorithm}[th]
\caption{$\textsc{Affine-Sensitivity}(P,w,p, k,\eps)$; see Theorem~\ref{lemma:pca}}\label{alg:tight}
{\begin{minipage}{\textwidth}
\begin{tabbing}
\textbf{Input:} \quad\=A weighted set $(P,w)$ of $n$ points in $\mathbb{R}^d$, a point $p\in P$,  \\\> an integer $k\in\br{0,\cdots,d-1}$, and an error parameter $\eps\in(0,1)$. \\
\textbf{Output:} An additive $\eps$-approximation $\tilde{s}$ to the sensitivity $s(p)$ of $p$.
\end{tabbing}\end{minipage}}

\nl$\psi := \left(\frac{\eps}{2}\right)^2$

\nl$\displaystyle{r := 1+ \max_{p\in P}\frac{D(p,0_d)}{\psi\eps^2}}$\label{lineneeded}

\nl$p' := (p\mid r)$  \quad\quad\tcp{$p'\in R^{d+1}$ is a concatenation of $p\in\ensuremath{\mathbb{R}}^d$ and $r\in\ensuremath{\mathbb{R}}$.}

\For {every $q \in P$}{

	 $q':=(q\mid r)$

	 $w(q') := w(q)$

}

$P' := \br{ q' \mid q \in P }$\label{definp}

$s' := \textsc{Non-Affine-Sensitivity}(P',w,p', k+1,\eps)$ \quad\quad\tcp{See Algorithm~\ref{alg:tight1}.}\label{defs}

$\tilde{s}:= s'+80\eps $

\Return $\tilde{s}$
\end{algorithm}

\begin{theorem} \label{lemma:pca} Let $(P,w)$ be a weighted set of $n$ points in $\ensuremath{\mathbb{R}}^d$, $p\in P$, $\eps\in(0,\frac{1}{12}]$ be an error parameter, and let $k\in \br{0,\cdots,d-1}$ be an integer. Let $\mathcal{S}^{A}_{d}$ denote the set of all affine $k$-subspaces in $\ensuremath{\mathbb{R}}^d$, $\displaystyle{s(p)=\sup_{S\in \mathcal{S}^{A}_{d}}\frac{w(p) D(p,S)}{\displaystyle{\sum_{q\in P}} w(q) D(q,S)}}$ denote the sensitivity of $p$ in the query space $(P,w,\mathcal{S}^{A}_{d},D)$, and let $\tilde{s}$ be the output of a call to $\textsc{Affine-Sensitivity}(P,w,p, k,\eps)$; See Algorithm~\ref{alg:tight}. Then
\[
s(p)\leq \tilde{s}\leq s(p)+161\eps.
\]
\end{theorem}

\begin{proof}
Let $S\in \mathcal{S}^{A}_{d}$. Notice that in this theorem $z=2$ and $D(p,S)=D_2(p,S)=\mathrm{dist}^2(p,S)$. Let  $\mathcal{S}_{d+1}$ denote the set of all non-affine $(k+1)$-subspaces of $\ensuremath{\mathbb{R}}^{d+1}$ and let $\displaystyle{s(p') = \sup_{S'\in \mathcal{S}_{d+1}}\frac{w(p)D(p',S')}{\displaystyle{\sum_{q\in P}} w(q) D(q',S')}}$. First by the definition of the set $P'$ at Line~\ref{definp} of Algorithm~\ref{alg:tight} and using Lemma~\ref{affinetononaffinecor} we have that for every point $p\in P$ and its corresponding $p'=(p\mid r)\in P'$ the following hold
\begin{align*}
 |s(p') - s(p)| \leq 16(2^z +1) \eps \cdot s(p') \leq   80\eps \cdot s(p') \leq 80\eps,
\end{align*}
where the last inequality holds since the sensitivity is always bounded by $1$ (i.e., $s(p')\leq1$). From the previous inequality we get
\begin{align}
s(p)-80 \eps \leq s(p') \leq s(p) + 80\eps. \label{pahes1}
\end{align}
 Let $s'$ be the output of a call to $\textsc{Non-Affine-Sensitivity}(P',w,p', k+1,\eps)$ as defined at Line~\ref{defs} of Algorithm~\ref{alg:tight}. By Lemma~\ref{lemma:svd} we have
\begin{align}
s(p') \leq s' \leq s(p') + \eps. \label{pahes2}
\end{align}
Combining~\eqref{pahes1} and~\eqref{pahes2} yields
$$s(p) \leq s(p')+80\eps \leq  s' + 80\eps \leq s(p') + 81\eps \leq s(p) + 161\eps, $$
where the first inequality holds by adding $80\eps$ to both sides of the left hand side inequality in~\eqref{pahes1}, the second and the third inequalities are by~\eqref{pahes2}, and the third holds by adding $81\eps$ to both side of the right hand side inequality in~\eqref{pahes1}. Considering the returned value $\tilde{s} = s' + 80\eps$  proves the theorem as,
$$s(p) \leq \tilde{s}\leq s(p) + 161\eps.$$
\end{proof}

\section{Experimental Results}\label{emp}
In this section we run benchmarks on real-world databases and compare our sampling algorithm with existing ones.

\paragraph{Algorithms. } We implemented the following sampling algorithms (distributions) for computing a coreset of $n$ points where every point is sampled with probability that is proportional to: (i) $1/n$ (uniform), (ii) existing sensitivity upper bound which is the sensitivity sampling algorithm of~\cite{DBLP:journals/corr/abs-1209-4893} that is mentioned at ``Introduction" section, and (iii) our tight sensitivity bound (Algorithm~\ref{alg:tight1} for SVD, and Algorithm~\ref{alg:tight} for PCA).

\paragraph{Software and Hardware. }
 We implemented those algorithms in Python 3.6 using the libraries Numpy~\cite{oliphant2006guide} and Scipy~\cite{SciPy}. We then run experimental results that we summarize in this section.  The experments where done using Intel i7-6850K CPU @ 3.60GHZ and 64GB RAM.

\subsection{Experimental Results for $k$-PCA}\label{pcaexp}
\paragraph{Datasets. } We used the following two datasets from~\cite{anguita2013public}:
(i) Gyroscope data, which we call  ``Gyro" in our graphs.
(ii) Embedded accelerometer data ($3$-axial linear accelerations) which we call ``Acc" in our graphs.
The data sets are resulted from experiments that have been carried out with a group of $30$ volunteers within an age bracket of $19$-$48$ years. Each person performed six activities (WALKING, WALKING UPSTAIRS, WALKING DOWNSTAIRS, SITTING, STANDING, LAYING) while wearing a smartphone (Samsung Galaxy S II) on the waist. Using its embedded gyroscope, $3$-axial angular velocities were recorded, at a constant frequency of $50$Hz. The experiments have been video-recorded to label the data manually.
Data was collected from $n=7352$ measurements. Each instance consists of measurements from $3$ dimensions, $x$, $y$, $z$, each in a size of $d=128$.  The results are those the corresponding $3$ datasets.

\paragraph{The experiment. } We ran Algorithms (i)-(iii) on the above datasets in order to compute sensitivities and sample coresets of variant sizes between $1000$ to $7000$. For each coreset, we computed the sum $OPT_k(A)$ of squared distances from the rows of the input matrix $A$ to the affine $k$-subspace that minimizes this sum. We then computed this sum to the optimal solution $OPT_k(C)$ on the coreset (to the rows of $A$). The approximation error $\varepsilon\in (0,1)$ was then defined to be $1-OPT_k(C)/OPT_k(A)$.

We used two values of $k= 5$ and $k=10$, and run each experiment 50 times.
 Results for the gyroscope data  (dataset (i)) are presented in Fig.~\ref{F1} and results for the accelerometer data (dataset (ii)) are presented in Fig.~\ref{F2}.

 	\begin{figure}[h!]
		\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/gyroxk5.png}

	\end{subfigure}
		\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/gyroyk5.png}

	\end{subfigure}
	\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/gyroxk5.png}

	\end{subfigure}

	\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/gyroxk10.png}

	\end{subfigure}
		\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/gyroyk10.png}

	\end{subfigure}
	\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/gyrozk10.png}

	\end{subfigure}

\caption{\small\label{F1} Experimental results of Subsection~\ref{pcaexp} for the gyroscope data.}
\end{figure}

	\begin{figure}[h!]
		\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/accxk5.png}

	\end{subfigure}
		\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/accyk5.png}

	\end{subfigure}
	\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/acczk5.png}

	\end{subfigure}

	\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/accxk10.png}

	\end{subfigure}
		\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/accyk10.png}

	\end{subfigure}
	\begin{subfigure}[h]{0.31\textwidth}
		\centering
		\includegraphics[scale=0.31]{small_data/acczk10.png}

	\end{subfigure}

\caption{\small\label{F2}Experimental results of Subsection~\ref{pcaexp} for the accelerometer data.}
\end{figure}

\subsection{Experimental Results for $k$-SVD}\label{svdexp}

\paragraph{Dataset}
We downloaded the document-term matrix of the English Wikipedia from~\cite{wic2019}, a sparse matrix of $4,624,611$ rows that correspond to documents, and $100$k columns (the dictionary of the $100$k most common words in Wikipedia~\cite{dic2012}). The entry in the $i$th row and $j$th column of this matrix is the number of how many appearances word number $j$ has in article number $i$. We call this data set ``Wiki" in our graphs.

\paragraph{Handling large data.} To handle this large dataset in memory, we maintain the coreset for the streaming set of rows, one by one, via the common merge and reduce tree that is usually used for this purpose; see e.g.\cite{feldman2010coresets} for details.

\paragraph{The experiment. }
We ran Algorithms (i)-(iii) on the document-term matrix of the English Wikipedia dataset in order to compute sensitivities and sample coresets of variant sizes between $1000$ to $7000$. For each coreset, we computed the sum $OPT_k(A)$ of squared distances from the rows of the input matrix $A$ to the (non-affine) $k$-subspace that minimizes this sum. We then computed this sum to the optimal solution $OPT_k(C)$ on the coreset (to the rows of $A$). The approximation error $\varepsilon\in (0,1)$ was then defined to be $1-OPT_k(C)/OPT_k(A)$.
 We ran with different values of $k$: $1$,$3$,$9$ and $11$.

\begin{figure}[h!]
	\begin{subfigure}[h]{0.5\textwidth}
		\centering
		\includegraphics[scale=0.5]{big_data_wiki/k=1_wiki1.png}

	\end{subfigure}
	\begin{subfigure}[h]{0.5\textwidth}
		\centering
		\includegraphics[scale=0.5]{big_data_wiki/k=3_wiki1.png}

	\end{subfigure}
	\begin{subfigure}[h]{0.5\textwidth}
		\centering
		\includegraphics[scale=0.5]{big_data_wiki/k=9_wiki1.png}

	\end{subfigure}
	\begin{subfigure}[h]{0.5\textwidth}
		\centering
		\includegraphics[scale=0.5]{big_data_wiki/k=11_wiki1.png}

	\end{subfigure}

\caption{\small\label{F5}Experimental results of Subsection~\ref{svdexp} for the document-term matrix of the English Wikipedia dataset.}
\end{figure}

\subsection{Conclusions and open problems}
We presented algorithms to compute exact sensitivity bounds for the $k$-SVD query spaces. Since the size of the coresets depends on the total sensitivity, we obtained coresets of size smaller and data dependent compared to existing worst-case upper bounds. We then suggested a generic reduction that enables us to generate tight sensivities also for the $k$-PCA problem (for affine $k$-subspaces).

Our experimental results show that our coreset indeed always smaller in practice.
We hope that the presented approach and open code would help to compute tight sensitivities for many other problems such as $k$-clustering, and other machine/deep learning problems.

\clearpage

\clearpage
\bibliographystyle{alpha}


\section{Introduction}
In the Standard Model (SM) it is assumed that the only difference between the three generations of leptons is their mass, and that their gauge couplings are the
same. This assumption, called lepton
universality (LU), can be tested by measuring the ratio of decay rates, which ensures that many experimental and theoretical uncertainties are cancelled in the ratio.
One type of these LU measurements is performed using semileptonic {\ensuremath{\PB}}\xspace decays of the form
${\ensuremath{\Pb}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\Pc}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$, commonly known as measurements of {\ensuremath{\mathcal{R}(\Hc)}}\xspace. This is defined as
\begin{equation}
    {\ensuremath{\mathcal{R}(\Hc)}}\xspace = \frac{\BF({\ensuremath{H_b}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{H_c}}\xspace{\ensuremath{\Ptau^-}}\xspace{\ensuremath{\neub_\tau}}\xspace)}{\BF({\ensuremath{H_b}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{H_c}}\xspace\ell^-{\ensuremath{\neub_\ell}}\xspace)} \, ,
\end{equation}
where {\ensuremath{H_b}}\xspace and {\ensuremath{H_c}}\xspace are a {\ensuremath{\Pb}}\xspace and {\ensuremath{\Pc}}\xspace hadron, respectively, and $\ell$ is either
an electron or muon.

Measurements of {\ensuremath{\mathcal{R}(\Hc)}}\xspace have been performed by the \mbox{LHCb}\xspace, \mbox{Belle}\xspace and \mbox{BaBar}\xspace experiments.
For {\ensuremath{\mathcal{R}(\D)}}\xspace, the average of the measured value of {\ensuremath{\mathcal{R}(\D)}}\xspace is
$0.349\pm 0.027\ensuremath{\mathrm{\,(stat)}}\xspace \pm 0.015\ensuremath{\mathrm{\,(syst)}}\xspace$ \cite{Lees:2012xj,Huschle:2015rga,Abdesselam:2019dgh}.
The predicted value for {\ensuremath{\mathcal{R}(\D)}}\xspace, assuming isospin symmetry, is
${\ensuremath{\mathcal{R}(\D)}}\xspace = {\ensuremath{\mathcal{R}(\Dp)}}\xspace = {\ensuremath{\mathcal{R}(\Dz)}}\xspace = 0.299 \pm 0.003$~\cite{Bigi:2016mdz,Bernlochner:2017jka,Jaiswal:2017rve,Aoki:2019cca}.
Even though {\ensuremath{\mathcal{R}(\D)}}\xspace differs from the SM prediction by only 1.4$\sigma$, the deviation from the SM of the
combined {\ensuremath{\mathcal{R}(\D)}}\xspace and {\ensuremath{\mathcal{R}(\Dstar)}}\xspace is about 3.1$\sigma$~\cite{HFLAV16}.

Radiative corrections were long thought to be negligible
at the level of precision of measurements and predictions of {\ensuremath{\mathcal{R}(\D)}}\xspace.
Recently, however, de Boer et al.~\cite{deBoer:2018ipi} presented a new evaluation of the
long-distance electromagnetic (QED) contributions
to ${\ensuremath{\Bbar{}^0}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^+}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$ and ${\ensuremath{\Bub}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^0}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$ decays, where ${\ensuremath{\ell^-}}\xspace = \mun, {\ensuremath{\Ptau^-}}\xspace$.
These corrections are different for \mun and {\ensuremath{\Ptau^-}}\xspace decays, such that they do not cancel in {\ensuremath{\mathcal{R}(\D)}}\xspace.
A proper evaluation of the radiative corrections can alter SM predictions and increase their uncertainty.

In measurements of {\ensuremath{\mathcal{R}(\D)}}\xspace, radiative corrections are simulated using the \mbox{\textsc{Photos}}\xspace package~\cite{Golonka:2005pn,Golonka:1379813}.
These proceedings, which are a summary of the studies described in Ref.~\cite{Cali:2019nwp}, show the difference between the QED
corrections simulated in \mbox{\textsc{Photos}}\xspace and those predicted by Ref.~\cite{deBoer:2018ipi}.
They describe a study on the effects of under- or overestimating these corrections in simulation on measurements of {\ensuremath{\mathcal{R}(\D)}}\xspace.

\section{Radiative corrections in P{\small HOTOS}}
\label{sec:comparePhotos}

\mbox{\textsc{Photos}}\xspace~\cite{Golonka:2005pn,Golonka:1379813} is a universal MC algorithm that simulates effects of QED corrections.
The corrections simulated by \mbox{\textsc{Photos}}\xspace have successfully been tested for {\ensuremath{\PW}}\xspace, {\ensuremath{\PZ}}\xspace, and {\ensuremath{\PB}}\xspace decays
and should be tested for every type of measurements individually, especially when high precision is needed.
Unlike Ref.~\cite{deBoer:2018ipi}, \mbox{\textsc{Photos}}\xspace does not include Coulomb corrections. These corrections concern the enhancement of decay rates
due to the interaction of two charged particles and are therefore relevant for the {\ensuremath{\D^+}}\xspace decay, but not for the {\ensuremath{\D^0}}\xspace decay.

Recent versions of \mbox{\textsc{Photos}}\xspace include multi-photon emissions as well as interference between final-state photons, whereas Ref.~\cite{deBoer:2018ipi}
also includes the interference between initial- and final state photons.
The calculation by de Boer et al. in Ref.~\cite{deBoer:2018ipi} is valid in the regime
in which the maximum energy of the radiated photons is smaller than the lepton mass, the muon mass in this case. \mbox{\textsc{Photos}}\xspace also includes photon
emission with higher energies. Neither Ref.~\cite{deBoer:2018ipi} nor \mbox{\textsc{Photos}}\xspace include the emission of photons depending on the hadronic structure.
These so called structure-dependent photons impact the spin of the decay particle and may also interfere with bremsstrahlung photons.

\begin{figure*}[t]
    \includegraphics[width=0.32\linewidth]{figure1a.pdf}
    \includegraphics[width=0.32\linewidth]{figure1b.pdf}
    \includegraphics[width=0.32\linewidth]{figure1c.pdf}
  \vspace{-0.2cm}
  \caption{
  Radiative corrections to the branching ratios of
  $\bar{B}^0\ensuremath{\rightarrow}\xspace D^+\ell^-\bar{\nu_\ell}$ (left) and $B^-\ensuremath{\rightarrow}\xspace D^0\ell^-\bar{\nu_\ell}$ (middle) decays, as a function of {\ensuremath{E_{\rm{max}}}}\xspace.
  The long-distance QED corrections to {\ensuremath{\mathcal{R}(\Dp)}}\xspace (orange) and {\ensuremath{\mathcal{R}(\Dz)}}\xspace
  (violet) as a function of {\ensuremath{E_{\rm{max}}}}\xspace (right).
  The plots are obtained from simulated data (solid
  lines) and from Ref.~\cite{deBoer:2018ipi} (dashed lines).}
  \label{fig:RD0plots}
\end{figure*}

\begin{figure*}
    \includegraphics[width=0.32\linewidth]{figure2a.pdf}
    \includegraphics[width=0.32\linewidth]{figure2b.pdf}
    \includegraphics[width=0.32\linewidth]{figure2c.pdf}
 \vspace{-0.2cm}
  \caption{
  Radiative corrections to the branching ratios of
  $\bar{B}^0\ensuremath{\rightarrow}\xspace D^+\ell^-\bar{\nu_\ell}$ (left) and {\ensuremath{\mathcal{R}(\Dp)}}\xspace (middle)
  in the case that no Coulomb correction is applied. The plots are
  obtained from simulated data (solid lines) and from
  Ref.~\cite{deBoer:2018ipi} (dashed lines). The plot on the right shows the ratio
  ${\ensuremath{\delta_{\rm{QED}}}}\xspace({\ensuremath{\mathcal{R}(\Dp)}}\xspace)/{\ensuremath{\delta_{\rm{QED}}}}\xspace({\ensuremath{\mathcal{R}(\Dz)}}\xspace)$.}
  \label{fig:noCoulomb}
\end{figure*}

To study the differences between Ref.~\cite{deBoer:2018ipi} and \mbox{\textsc{Photos}}\xspace, four samples (${\ensuremath{\Bbar{}^0}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^+}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$ and ${\ensuremath{\Bub}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^0}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$,
where ${\ensuremath{\ell^-}}\xspace = \mun, {\ensuremath{\Ptau^-}}\xspace$) with three million {\ensuremath{\PB}}\xspace meson decays are generated by \mbox{\textsc{Pythia}}\xspace 8~\cite{Sjostrand:2006za,Sjostrand:2007gs}.
The {\ensuremath{\PB}}\xspace mesons decays are simulated by \mbox{\textsc{EvtGen}}\xspace~\cite{Lange:2001uf}, including the QED corrections by \mbox{\textsc{Photos}}\xspace v.3.56, with the ``option with interference'' switched on.
For both the ${\ensuremath{\PB}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\PD}}\xspace\mu{\bar{\nu}}_\mu$ and ${\ensuremath{\PB}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\PD}}\xspace\tau{\bar{\nu}}_\tau$ decays considered, the \verb+HQET2+ model is used with parameters from Ref.~\cite{HFLAV16}.

The four-momentum of the total radiated photons, $\mbox{$p$}\xspace_{\gamma}$, is defined as
\begin{equation}
    \mbox{$p$}\xspace_{\gamma} = \mbox{$p$}\xspace_{B} - \left(\mbox{$p$}\xspace_{{\ensuremath{\PD}}\xspace} + \mbox{$p$}\xspace_{{\ensuremath{\ell^-}}\xspace} + \mbox{$p$}\xspace_{{\ensuremath{\neub_\ell}}\xspace} \right) \, ,
\end{equation}
where $\mbox{$p$}\xspace_{B}$, $\mbox{$p$}\xspace_{{\ensuremath{\PD}}\xspace}$, $\mbox{$p$}\xspace_{{\ensuremath{\neub_\ell}}\xspace}$, and $\mbox{$p$}\xspace_{{\ensuremath{\ell^-}}\xspace}$ are the four-momenta of the
{\ensuremath{\PB}}\xspace, {\ensuremath{\PD}}\xspace, {\ensuremath{\ell^-}}\xspace and {\ensuremath{\neub_\ell}}\xspace, respectively.
In agreement with Ref.~\cite{deBoer:2018ipi}, the radiation of the {\ensuremath{\PD}}\xspace decay products is not taken into account.
The total energy of the radiated photons, $E_{\gamma}$, is computed in the {\ensuremath{\PB}}\xspace rest frame.

The variable {\ensuremath{E_{\rm{max}}}}\xspace is defined as the maximum value that $E_{\gamma}$ is allowed to have to consider ${\ensuremath{\PB}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\PD}}\xspace\ell{\bar{\nu}}_\ell(\gamma)$ signal rather than background.
The QED correction, {\ensuremath{\delta_{\rm{QED}}}}\xspace, is given by the relative variation of the branching ratio when events with total radiated energy greater than {\ensuremath{E_{\rm{max}}}}\xspace are discarded, calculated as
\begin{equation}
{\ensuremath{\delta_{\rm{QED}}}}\xspace = \frac{\int_{0}^{{\ensuremath{E_{\rm{max}}}}\xspace} N(E_{\gamma} ) dE_{\gamma} }{\int_{0}^{\infty} N(E_{\gamma}) dE_{\gamma}} - 1\, .
\end{equation}
Here, $N(E_{\gamma} )$ is the distribution of events with $E_{\gamma}$.
The considered energy range is up to 100~MeV, covering the majority of radiative photons generated by \mbox{\textsc{Photos}}\xspace.

Figure~\ref{fig:RD0plots} shows comparisons between radiative corrections from \mbox{\textsc{Photos}}\xspace
and Ref.~\cite{deBoer:2018ipi} to the ${\ensuremath{\Bbar{}^0}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^+}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$ (left) and ${\ensuremath{\Bub}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^0}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$ (middle) branching ratios.
These plots show differences of up to 2\% for the {\ensuremath{\Bbar{}^0}}\xspace decays, and $0.5-1\%$ for {\ensuremath{\Bub}}\xspace decays.
This effect does not cancel even in the ratios of branching fractions. This is clearly visible in Fig.~\ref{fig:RD0plots} (right),
where radiative corrections on {\ensuremath{\mathcal{R}(\D)}}\xspace, $\delta_{QED}(\mathcal{R})$, are shown as a function of {\ensuremath{E_{\rm{max}}}}\xspace.
\mbox{\textsc{Photos}}\xspace predicts a QED correction of 0.5\% lower than Ref.~\cite{deBoer:2018ipi} in {\ensuremath{\mathcal{R}(\Dp)}}\xspace, and 0.5\% higher in {\ensuremath{\mathcal{R}(\Dz)}}\xspace.

A significant part of the radiative corrections in Ref.~\cite{deBoer:2018ipi} originates
from Coulomb interactions, which are not included in \mbox{\textsc{Photos}}\xspace.
Light leptons typically have a Coulomb correction of 1.023~\cite{Atwood:1989em},
whereas the {\ensuremath{\Ptau^-}}\xspace leptons in the ${\ensuremath{\Bbar{}^0}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^+}}\xspace{\ensuremath{\Ptau^-}}\xspace{\ensuremath{\neub_\tau}}\xspace$ decay have Coulomb corrections between 2.5\% and 5.0\%.
The QED corrections from \mbox{\textsc{Photos}}\xspace for the {\ensuremath{\D^+}}\xspace decay mode
are compared with predictions not including the Coulomb correction as provided by Ref.~\cite{deBoer:2018ipi}.
This reduces the difference of the corrections to the branching ratios between \mbox{\textsc{Photos}}\xspace and the theoretical calculations
to about 1\% and brings the corrections on {\ensuremath{\mathcal{R}(\Dp)}}\xspace in close agreement, as is shown in Fig.~\ref{fig:noCoulomb}.

Figure~\ref{fig:noCoulomb} (right) shows the ratio of QED corrections on {\ensuremath{\mathcal{R}(\Dp)}}\xspace over those on {\ensuremath{\mathcal{R}(\Dz)}}\xspace. Both \mbox{\textsc{Photos}}\xspace and the calculation in Ref.~\cite{deBoer:2018ipi} without Coulomb corrections conserve isospin symmetry ({\ensuremath{\delta_{\rm{QED}}}}\xspace values for {\ensuremath{\mathcal{R}(\Dp)}}\xspace and {\ensuremath{\mathcal{R}(\Dz)}}\xspace agree within the errors), while Coulomb corrections introduce an isospin-breaking term.

\section{Effects on LHCb-like analysis}
\label{sec:dummyAnalysis}

The comparison between Ref.~\cite{deBoer:2018ipi} and \mbox{\textsc{Photos}}\xspace can be made only for soft photons with energies up to 100~MeV. For higher energies, no
calculations are available. However, \mbox{\textsc{Photos}}\xspace generates also photons with higher energies in ranges where structure-dependent photons are relevant.
These are used to study the effects of under- or overestimating radiative corrections in simulation for a measurement of {\ensuremath{\mathcal{R}(\D)}}\xspace in an \mbox{LHCb}\xspace-like environment.

A study is performed with the same data sets as described in the previous section by making a template fit to the
${\ensuremath{\PB}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\PD}}\xspace\mun{\ensuremath{\neub_\mu}}\xspace$ and ${\ensuremath{\PB}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\PD}}\xspace{\ensuremath{\Ptau^-}}\xspace{\ensuremath{\neub_\tau}}\xspace$ components. This fit uses the same fit variables as \mbox{LHCb}\xspace's muonic {\ensuremath{\mathcal{R}(\Dstar)}}\xspace analysis~\cite{Aaij:2015yra}.
These are the muon energy computed in the {\ensuremath{\PB}}\xspace meson rest frame, {\ensuremath{E_{\mu}}}\xspace; the missing mass squared, ${\ensuremath{m_{\rm miss}^2}}\xspace = (\mbox{$p$}\xspace_{{\ensuremath{\PB}}\xspace}-\mbox{$p$}\xspace_{{\ensuremath{\PD}}\xspace}-\mbox{$p$}\xspace_{\mu})^2$; and the squared four-momentum transferred to the lepton system, ${\ensuremath{q^2}}\xspace=(\mbox{$p$}\xspace_{{\ensuremath{\PB}}\xspace}-\mbox{$p$}\xspace_{{\ensuremath{\PD}}\xspace})^2$. Simulated data samples are created from a mixture of ${\ensuremath{\PB}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\PD}}\xspace\mun{\ensuremath{\neub_\mu}}\xspace$ and ${\ensuremath{\PB}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\PD}}\xspace{\ensuremath{\Ptau^-}}\xspace{\ensuremath{\neub_\tau}}\xspace$ decays, with radiative corrections generated by \mbox{\textsc{Photos}}\xspace. Here {\ensuremath{\mathcal{R}(\D)}}\xspace is assumed to be 0.3 as predicted by the SM. No backgrounds
are considered.

The fits are performed with templates that are created under the hypothesis that there is no radiation $E_\gamma$ above {\ensuremath{E_{\rm{max}}}}\xspace. Five values of {\ensuremath{E_{\rm{max}}}}\xspace,
100, 300, 500, 800 and 1500~MeV, are chosen for this study.
Fitting the templates to the simulated data sample with no cuts on radiation yields an estimate of the possible bias on {\ensuremath{\mathcal{R}(\D)}}\xspace.
This indicates the importance of simulating $E_\gamma$ in the high-energy region. Note that \mbox{LHCb}\xspace analyses do not cut on radiative energy explicitly, but that
implicitly applied cuts could alter this bias.

 \begin{figure}[t]
  \includegraphics[width=0.85\linewidth]{figure3.pdf}
  \vspace{-0.4cm}
  \caption{
   Ratios of the cut {\ensuremath{m_{\rm miss}^2}}\xspace distribution over the default {\ensuremath{m_{\rm miss}^2}}\xspace distribution for the  $\bar{B}^0\ensuremath{\rightarrow}\xspace D^+\ell^-\bar{\nu_\ell}$
  decays, for the various cuts on {\ensuremath{E_{\rm{max}}}}\xspace for the \mun decay mode.}
  \label{fig:mmiss_emaxcuts}
\end{figure}

The acceptance of the \mbox{LHCb}\xspace detector is mimicked by applying selection requirements following Ref.~\cite{Ciezarek:2016lqu}.
The production and decay vertices are smeared to simulate the detector resolution and a cut on the flight distance is applied
to resemble the trigger selection. At \mbox{LHCb}\xspace, the {\ensuremath{\PB}}\xspace meson momentum cannot be reconstructed due to the missing neutrino.
Therefore, as in Ref.~\cite{Aaij:2015yra}, the momentum of the {\ensuremath{\PB}}\xspace in the $z$ direction, $(\mbox{$p$}\xspace_{{\ensuremath{\PB}}\xspace})_z$, is approximated as
$(\mbox{$p$}\xspace_{{\ensuremath{\PB}}\xspace})_z=(m_{{\ensuremath{\PB}}\xspace}/m_{\rm vis})(\mbox{$p$}\xspace_{\rm vis})_z$, where $m_{{\ensuremath{\PB}}\xspace}$ is the {\ensuremath{\PB}}\xspace mass,
and $m_{\rm vis}$ and $(\mbox{$p$}\xspace_{\rm vis})_z$ are the mass and momentum in the
beam direction of the visible decay products of the {\ensuremath{\PB}}\xspace meson, respectively.

The simulated sample size is based on an estimate of the yields for the Run II data-taking period using
the reconstruction efficiencies from Ref.~\cite{Aaij:2015yra}, the {\ensuremath{\PB}}\xspace production cross-section at 13~TeV,
and branching fractions. This results in yields of data samples of $1.0\times10^6$ and $0.5\times10^5$ for the
${\ensuremath{\Bbar{}^0}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^+}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$ decays, and $4.4\times 10^5$ and $2.3\times10^4$ for
the ${\ensuremath{\Bub}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^0}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$ decays, where the first yield represents the \mun sample, and the second the {\ensuremath{\Ptau^-}}\xspace sample.

The values of {\ensuremath{\mathcal{R}(\D)}}\xspace are determined from the fitted yields as well as the reconstruction efficiencies
$\varepsilon_{\mu}$ and $\varepsilon_{\tau}$ for the \mun and {\ensuremath{\Ptau^-}}\xspace samples using
\begin{equation}
    {\ensuremath{\mathcal{R}(\D)}}\xspace = \frac{f_{\tau}}{1-f_{\tau}} \frac{\varepsilon_{\mu}}{\varepsilon_{\tau}} \, .
\end{equation}
It is found that for this specific case study, the ratio of efficiencies is not affected by the cuts on $E_{\gamma}$.
Combining the efficiencies with the fitted yields, the resulting values of {\ensuremath{\mathcal{R}(\Dp)}}\xspace as a function of the cut on {\ensuremath{E_{\rm{max}}}}\xspace are shown in Fig.~\ref{fig:RD_Emax},
which shows a dependence on {\ensuremath{E_{\rm{max}}}}\xspace. From here it is clear that there is a significant effect in
underestimating the QED radiative corrections which could be up to 0.02 for both {\ensuremath{\mathcal{R}(\Dp)}}\xspace and {\ensuremath{\mathcal{R}(\Dz)}}\xspace values, corresponding to a relative bias of 7.5\%.
This is due to the change in template shapes when applying cuts on {\ensuremath{E_{\rm{max}}}}\xspace which is most clearly visible in the {\ensuremath{m_{\rm miss}^2}}\xspace distribution of the semimuonic decay, shown in
Fig.~\ref{fig:mmiss_emaxcuts}. Since the template shapes of the semitauonic decays do not change significantly, this effect does not cancel in the ratio {\ensuremath{\mathcal{R}(\Dp)}}\xspace.

\begin{figure}[h]
  \includegraphics[width=0.85\linewidth]{figure4.pdf}
  \vspace{-0.4cm}
  \caption{Measured values of {\ensuremath{\mathcal{R}(\Dp)}}\xspace in a simplified \mbox{LHCb}\xspace-like analysis, as a function of {\ensuremath{E_{\rm{max}}}}\xspace.
  The blue band shows the fit result obtained with the same templates used to generate the pseudo-experiments.
  }
  \label{fig:RD_Emax}
\end{figure}

The results for {\ensuremath{\mathcal{R}(\Dz)}}\xspace analysis look very similar to those for {\ensuremath{\mathcal{R}(\Dp)}}\xspace and are therefore omitted from these proceedings.
In actual analyses, radiative corrections are present in data, and, at least partially, in simulation. Therefore this
analysis is likely to be an overestimate on the bias.

Also the Coulomb correction impacts the shape of the fit templates and thus the experimental results of {\ensuremath{\mathcal{R}(\Dp)}}\xspace.
This bias is evaluated by weighting each event in the ${\ensuremath{\Bbar{}^0}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^+}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$ decay by the Coulomb correction {\ensuremath{\Omega_{\rm{C}}}}\xspace.
Changes in the {\ensuremath{q^2}}\xspace, {\ensuremath{m_{\rm miss}^2}}\xspace and {\ensuremath{E_{\mu}}}\xspace distributions are shown in Fig.~\ref{fig:qsquare_Coulomb}. {\ensuremath{\Omega_{\rm{C}}}}\xspace is mostly constant for the \mun mode,
but for the {\ensuremath{\Ptau^-}}\xspace mode, a dependence on each of the three variables is shown. The above analysis is repeated while applying Coulomb corrections
to the simulated data sample and not on the fit templates, resulting in a relative shift on {\ensuremath{\mathcal{R}(\Dp)}}\xspace of about -1.0\%. No additional cuts on {\ensuremath{E_{\rm{max}}}}\xspace are applied.

\begin{figure*}[t]
    \includegraphics[width=0.32\linewidth]{figure5a.pdf}
    \includegraphics[width=0.32\linewidth]{figure5b.pdf}
    \includegraphics[width=0.32\linewidth]{figure5c.pdf}
  \vspace{-0.2cm}
  \caption{
  Coulomb corrections as a function of {\ensuremath{q^2}}\xspace, {\ensuremath{m_{\rm miss}^2}}\xspace, and {\ensuremath{E_{\mu}}}\xspace for the ${\ensuremath{\Bbar{}^0}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\D^+}}\xspace{\ensuremath{\ell^-}}\xspace{\ensuremath{\neub_\ell}}\xspace$ decays, where ${\ensuremath{\ell^-}}\xspace = \mun, {\ensuremath{\Ptau^-}}\xspace$.}
  \label{fig:qsquare_Coulomb}
\end{figure*}

\newpage
\section{Conclusions and recommendations}
\label{sec:Conclusions}

The QED corrections described in Ref.~\cite{deBoer:2018ipi} are not fully included in the \mbox{\textsc{Photos}}\xspace
package which is used to simulate these corrections in analyses of the \mbox{LHCb}\xspace, \mbox{Belle}\xspace, and \mbox{BaBar}\xspace experiments.
These different QED corrections affect the muonic and tauonic branching ratios
at the level of a few percent. While calculating the ratios {\ensuremath{\mathcal{R}(\D)}}\xspace, the differences largely cancel out
when Coulomb corrections are ignored. Coulomb corrections are not present in \mbox{\textsc{Photos}}\xspace and this results
in a discrepancy between Ref.~\cite{deBoer:2018ipi} and \mbox{\textsc{Photos}}\xspace of up to 1~\% on the ratio {\ensuremath{\mathcal{R}(\Dp)}}\xspace.

Coulomb interactions mainly affect the kinematics of tauonic decays, changing the shape of the
distributions used to determine their signal yields in an \mbox{LHCb}\xspace-like analysis. Not including these in simulated data can
result in a bias of around 1\% on measurements of {\ensuremath{\mathcal{R}(\D)}}\xspace.
Over- or underestimating radiative corrections can bias \mbox{LHCb}\xspace-like measurements up to 7\%, resulting
in a bias of 0.02 on {\ensuremath{\mathcal{R}(\Dp)}}\xspace.

These studies must be repeated for each analysis individually because they are
dependent on selection requirements as well as fit variables.
For \mbox{Belle}\xspace~II measurements, which have a better resolution on the kinematic variables~\cite{Kou:2018nap} than \mbox{LHCb}\xspace,
the effects could even be larger.
Additional calculations of QED corrections on for ${\ensuremath{\PB}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\PD}}\xspace\ell{\ensuremath{\neu_\ell}}\xspace$ decays, specifically those involving high-energy and
structure-dependent photons, are necessary in order to make measurements with higher precision.

\begin{acknowledgments}
These proceedings are a summary of the work presented in Ref.~\cite{Cali:2019nwp} and I am very grateful to my fellow authors
Stefano Cal\'{i}, Marcello Rotondo and Barbara Sciascia for the pleasant collaboration.
\end{acknowledgments}

\bigskip