Numerical identifiability and optimal experiment design

Let $p\ge 1$, let $\mathrm{\Theta }\subseteq {\mathbb{R}}^p$ be open, and let $\mathcal{X}$ be a finite nonempty set of experimental conditions. For each $x\in \mathcal{X}$, let $q_x\ge 1$, let $F_x:\mathrm{\Theta }\to {\mathbb{R}}^{q_x}$ be $C^1$ in the sense of Continuously differentiable maps, local inverses, and local diffeomorphisms, and let $W_x\in {\mathbb{R}}^{q_x\times q_x}$ be a symmetric positive-definite weight matrix, meaning $z^T{W}_{x}z>0$ for every nonzero $z\in {\mathbb{R}}^{q_x}$.

An approximate design is a function $\xi :\mathcal{X}\to [0,1]$ with ${\sum }_{x\in \mathcal{X}}{\xi }_x=1$, as in Approximate designs and linear-model optimality criteria. Its support is $\mathrm{supp}(\xi )=\{x\in \mathcal{X}:{\xi }_x>0\}$. Define the joint observation map by $F_{\xi }(\theta )=(F_x(\theta ){)}_{x\in \mathrm{supp}(\xi )},$ where the finite product is ordered by any fixed listing of $\mathrm{supp}(\xi )$. Its stacked sensitivity Jacobian at $\theta$ is the vertical concatenation $J_{\xi }(\theta )$ of the matrices $J_x(\theta )=D{F}_x(\theta )$ for $x\in \mathrm{supp}(\xi )$. The associated weighted information matrix is $M_{\theta }(\xi )={\sum }_{x\in \mathcal{X}}{\xi }_x{J}_x(\theta {)}^T{W}_x{J}_x(\theta )\mathrm{.}$ The matrix is symmetric positive semidefinite because each summand has nonnegative quadratic form.

Because the support is finite, $F_{\xi }$ is $C^1$ with derivative $J_{\xi }(\theta )$: the finite collection of component first-order remainders has a stacked Euclidean norm that is still $o(\mathrm{\parallel }h{\mathrm{\parallel }}_2)$, and the stacked derivative is continuous because each $D{F}_x$ is continuous.

For a fixed design $\xi$, the observation map is locally identifiable at ${\theta }_0\in \mathrm{\Theta }$ if it is injective on some open neighborhood of ${\theta }_0$ in $\mathrm{\Theta }$. It is globally identifiable on $\mathrm{\Theta }$ if it is injective on all of $\mathrm{\Theta }$. A matrix $J:{\mathbb{R}}^p\to {\mathbb{R}}^m$ has full column rank when $Jv=0$ implies $v=0$. Its smallest singular value is ${\sigma }_{\min }(J)={\inf }_{\mathrm{\parallel }v{\mathrm{\parallel }}_2=1}\mathrm{\parallel }Jv{\mathrm{\parallel }}_2$; thus ${\sigma }_{\min }(J)\ge \mu$ means $\mathrm{\parallel }Jv{\mathrm{\parallel }}_2\ge \mu \mathrm{\parallel }v{\mathrm{\parallel }}_2$ for every $v$.

The approximate design uses real weights. An exact $n$-run design instead restricts weights to multiples $n_x\mathrm{/}n$ with integer counts, and is a separate finite allocation problem in Approximate designs and linear-model optimality criteria.

Let $J,J_h\in {\mathbb{R}}^{m\times n}$, let $r=\min \{m,n\}$, and let ${\sigma }_1(A)\ge \cdots \ge {\sigma }_r(A)\ge 0$ denote the singular values of a real matrix $A$. Then for every $1\le i\le r$, $\mathrm{\mid }{\sigma }_i(J_h)-{\sigma }_i(J)\mathrm{\mid }\le \mathrm{\parallel }{J}_h-J{\mathrm{\parallel }}_2,$ where $\mathrm{\parallel }\cdot {\mathrm{\parallel }}_2$ is the Euclidean operator norm.

Let $p\ge 1$ and let $A,B\in {\mathbb{R}}^{p\times p}$ be symmetric. Put $E=B-A$ and $A_t=(1-t)A+tB$ for $0\le t\le 1$. Suppose there is $\alpha >0$ such that $v^T{A}_{t}v\ge \alpha \mathrm{\parallel }v{\mathrm{\parallel }}_2^2$ for every $v\in {\mathbb{R}}^p$ and every $t\in [0,1]$. Then $A_t$ is positive definite and $\mathrm{\mid }\log \det B-\log \det A\mathrm{\mid }\le \frac{p\mathrm{\parallel }E{\mathrm{\parallel }}_2}{\alpha }\mathrm{.}$

Let $p\ge 1$, let $\mathcal{X}$ be a finite nonempty candidate set, and let $f_x\in {\mathbb{R}}^p$ be the feature vector for each $x\in \mathcal{X}$. An approximate design $\xi$ and information matrix $M(\xi )={\sum }_x{\xi }_x{f}_x{f}_x^T$ use the convention of Approximate designs and linear-model optimality criteria. The same vectors may be scalar-observation sensitivities at a fixed parameter, as in Finite nonlinear observation designs, local identifiability, and weighted information, with a scalar precision absorbed into $f_x$. The procedure uses exact arithmetic.

Input: A finite nonempty $\mathcal{X}$, $p\ge 1$, its feature vectors $(f_x{)}_{x\in \mathcal{X}}$, an approximate initial design ${\xi }^{(0)}$, and a finite nonnegative integer budget $N$ counting weight updates.

Procedure: Form $M_0=M({\xi }^{(0)})$. If $M_0$ is not positive definite, return $({\xi }^{(0)},\text{undefined violation},\text{invalid initial design})$. Otherwise, for $k=0,1,\dots ,N$, with current ${\xi }^{(k)}$ and $M_k=M({\xi }^{(k)})$, compute $s_x^{(k)}=f_x^T{M}_k^{-1}{f}_x(x\in \mathcal{X}),v_k={\max }_{x\in \mathcal{X}}{s}_x^{(k)}-p\mathrm{.}$ If $v_k=0$, return $({\xi }^{(k)},v_k,\text{equivalence-certified})$. If $k=N$, return $({\xi }^{(k)},v_k,\text{budget-exhausted})$. Otherwise set ${\xi }_x^{(k+1)}={\xi }_x^{(k)}\frac{s_x^{(k)}}{p}(x\in \mathcal{X}),$ form $M_{k+1}={\sum }_x{\xi }_x^{(k+1)}{f}_x{f}_x^T$, and continue.

Finite-budget multiplicative design updates preserve the simplex and nonsingularity ↗ proves preservation of positive definiteness and the bounds of at most $N$ updates and $N+1$ sensitivity scans. A nonzero returned violation is reported as-is; it is not an exact D-optimality certificate. The exact integer allocation problem remains separate from this approximate-weight procedure.

Use the finite observation designs of Finite nonlinear observation designs, local identifiability, and weighted information. Suppose uniformly over $x\in \mathcal{X}$ that $\mathrm{\parallel }{J}_x{\mathrm{\parallel }}_2\le K$, $\mathrm{\parallel }{J}_{h,x}-J_x{\mathrm{\parallel }}_2\le \eta$, and $\mathrm{\parallel }{W}_x{\mathrm{\parallel }}_2\le w$, where $K,\eta ,w\ge 0$. For every approximate design $\xi$, $\mathrm{\parallel }{M}_h(\xi )-M(\xi ){\mathrm{\parallel }}_2\le \delta :=w(2K\eta +{\eta }^2),$ where $M_h(\xi )={\sum }_x{\xi }_x{J}_{h,x}^T{W}_x{J}_{h,x}$. If $M(\xi )$ satisfies $v^{T}M(\xi )v\ge \alpha \mathrm{\parallel }v{\mathrm{\parallel }}_2^2$ for every $v$ with some $\alpha >\delta$, then $v^T{M}_h(\xi )v\ge (\alpha -\delta )\mathrm{\parallel }v{\mathrm{\parallel }}_2^2$ for every $v$, so the computed information matrix remains positive definite.

For a $C^2$ residual map $r:U\to {\mathbb{R}}^m$ on an open $U\subseteq {\mathbb{R}}^n$, with $m,n\ge 1$, consider the unweighted nonlinear least-squares objective $\varphi (\theta )=\frac{1}{2}\mathrm{\parallel }r(\theta ){\mathrm{\parallel }}_2^2$. Then the matrix $J(\theta {)}^{T}J(\theta )$ is the Gauss–Newton term of the Hessian. The exact identity is ${\mathrm{\nabla }}^2\varphi (\theta )=J(\theta {)}^{T}J(\theta )+{\sum }_i{r}_i(\theta ){\mathrm{\nabla }}^2{r}_i(\theta )$ by The nonlinear least-squares gradient and Hessian formulas. Thus the two matrices agree at zero residual, but the information matrix defined in Finite nonlinear observation designs, local identifiability, and weighted information is not asserted to be the full Hessian at a nonzero-residual point.

Let $\mathrm{\Theta }\subseteq {\mathbb{R}}^p$ be open with $p\ge 1$, let $F_{\xi }:\mathrm{\Theta }\to {\mathbb{R}}^m$ be the joint observation map of Finite nonlinear observation designs, local identifiability, and weighted information, and fix ${\theta }_0\in \mathrm{\Theta }$. If its Jacobian $J_{\xi }({\theta }_0)$ has full column rank, then $F_{\xi }$ is injective on an open neighborhood of ${\theta }_0$; hence the model is locally identifiable at ${\theta }_0$ for this design.

For the finite nonlinear observation design of Finite nonlinear observation designs, local identifiability, and weighted information, fix $\theta \in \mathrm{\Theta }$ and an approximate design $\xi$. The matrix $M_{\theta }(\xi )$ is positive definite, meaning symmetric with $v^T{M}_{\theta }(\xi )v>0$ for every nonzero $v\in {\mathbb{R}}^p$, if and only if the stacked support sensitivity Jacobian $J_{\xi }(\theta )$ has trivial kernel, equivalently full column rank.

Let $F:\mathrm{\Theta }\to {\mathbb{R}}^m$ be $C^1$ on an open $\mathrm{\Theta }\subseteq {\mathbb{R}}^p$, let ${\theta }_0\in \mathrm{\Theta }$, and let $r_0>0$ be such that the closed Euclidean ball $\bar{B}({\theta }_0,r_0)$ lies in $\mathrm{\Theta }$. Suppose $\mathrm{\parallel }DF({\theta }_0)v{\mathrm{\parallel }}_2\ge \mu \mathrm{\parallel }v{\mathrm{\parallel }}_2$ for every $v\in {\mathbb{R}}^p$, where $\mu >0$, and suppose $DF$ is $L$-Lipschitz in operator norm on this ball, where $L\ge 0$. Define $r=r_0$ if $L=0$ and $r=\min \{r_0,\mu \mathrm{/}(2L)\}$ if $L>0$. If $y=F({\theta }_0)+e$ with $\mathrm{\parallel }e{\mathrm{\parallel }}_2\le \delta$, and $\widehat{\theta }\in \bar{B}({\theta }_0,r)$ satisfies $\mathrm{\parallel }F(\widehat{\theta })-y{\mathrm{\parallel }}_2\le \epsilon$, then $\mathrm{\parallel }\widehat{\theta }-{\theta }_0{\mathrm{\parallel }}_2\le \frac{2(\epsilon +\delta )}{\mu }\mathrm{.}$ The bound is conditional on the stated local ball and on the existence of such a fit; it does not assert that a fit exists.

For the algorithm Finite-budget multiplicative update for approximate D-optimal design, suppose the initial weights form an approximate design and its information matrix is positive definite. Then every update preserves nonnegative weights summing to one and positive definiteness of the information matrix. For every finite budget $N$, the procedure returns after at most $N$ updates and at most $N+1$ sensitivity scans.

Let $\mathcal{X}$ be a finite candidate set and let $\mathrm{\Xi }$ be a nonempty compact metric space of approximate designs on $\mathcal{X}$. For each $x$, let $J_x$ and $J_{h,x}$ be exact and computed sensitivities and let $W_x$ be symmetric positive definite. Suppose uniformly over $x\in \mathcal{X}$ that $\mathrm{\parallel }{J}_x{\mathrm{\parallel }}_2\le K$, $\mathrm{\parallel }{J}_{h,x}-J_x{\mathrm{\parallel }}_2\le \eta$, and $\mathrm{\parallel }{W}_x{\mathrm{\parallel }}_2\le w$, where $K,\eta ,w\ge 0$. Put $\delta =w(2K\eta +{\eta }^2)$. Suppose the criteria $\mathrm{\Phi }(\xi )=\log \det M(\xi )$ and ${\mathrm{\Phi }}_h(\xi )=\log \det M_h(\xi )$ are continuous on $\mathrm{\Xi }$, and for some $\alpha >0$ every matrix on every segment between $M(\xi )$ and $M_h(\xi )$ satisfies $v^{T}Av\ge \alpha \mathrm{\parallel }v{\mathrm{\parallel }}_2^2$. If ${\xi }^{\ast }$ maximizes $\mathrm{\Phi }$ and $\widehat{\xi }$ is a $\tau$-approximate maximizer of ${\mathrm{\Phi }}_h$, meaning ${\mathrm{\Phi }}_h(\widehat{\xi })\ge {\max }_{\xi \in \mathrm{\Xi }}{\mathrm{\Phi }}_h(\xi )-\tau$ with $\tau \ge 0$, then $\mathrm{\Phi }({\xi }^{\ast })-\mathrm{\Phi }(\widehat{\xi })\le \tau +\frac{2p\delta }{\alpha }\mathrm{.}$

Let $\mathcal{X}$ be a finite nonempty candidate set, let $p\ge 1$, and let $\xi$ be any returned nonsingular approximate design from Finite-budget multiplicative update for approximate D-optimal design. Put $M=M(\xi )$ and $s_x=f_x^T{M}^{-1}{f}_x,v={\max }_{x\in \mathcal{X}}{s}_x-p\mathrm{.}$ Then $v\ge 0$, every nonsingular approximate competitor $\eta$ satisfies $\log \det M(\eta )-\log \det M(\xi )\le v,$ and $v=0$ if and only if $\xi$ is D-optimal among nonsingular approximate designs.