Inversion-Free Evaluation of Nearest Neighbors in Method of Moments

TL;DR

This paper presents an inversion-free, scalable method for evaluating nearest neighbors in the method of moments, enhancing shape perturbation analysis and data mining applications.

Contribution

It extends topology sensitivity in method of moments by incorporating degrees-of-freedom reconstruction, enabling efficient parallelizable nearest neighbor evaluation.

Findings

Effective for small shape perturbation evaluation

Suitable for parallel computation

Useful for machine learning data mining

Abstract

A recently introduced technique of topology sensitivity in method of moments is extended by the possibility of adding degrees-of-freedom (reconstruct) into underlying structure. The algebraic formulation is inversion-free, suitable for parallelization and scales favorably with the number of unknowns. The reconstruction completes the nearest neighbors procedure for an evaluation of the smallest shape perturbation. The performance of the method is studied with a greedy search over a Hamming graph representing the structure in which initial positions are chosen from a random set. The method is shown to be effective data mining tool for machine learning-related applications.

Figures

Tables

Table 1. TABLE I: Comparison of shape optimization performance depending on the number of d-o-f . The table is extended by a column with the results for an 8 × 16 8 16 8\times 16 grid.

plate	$4\times 8$	$6\times 12$	$8\times 16$
d-o-f, $N$	$180$	$414$	$744$
runs, $I$	$5\cdot {10}^4$	$5\cdot {10}^4$	$1\cdot {10}^3$
comp. time, $T$ [s]	$2.4\cdot {10}^3$	$5.8\cdot {10}^4$	$1.2\cdot {10}^4$
evaluated shapes	$7.2\cdot {10}^8$	$3.9\cdot {10}^9$	$2.6\cdot {10}^8$
shapes per second	$3\cdot {10}^5$	$7\cdot {10}^4$	$2\cdot {10}^4$
comp. time per run, $T/I$ [s]	$4.8\cdot {10}^{-2}$	$1.2\cdot {10}^0$	$1.2\cdot {10}^1$
evaluated shapes per run	$1.4\cdot {10}^4$	$7.8\cdot {10}^4$	$2.6\cdot {10}^5$
$Q_{\min }/Q_{\mathrm{lb}}^{\mathrm{TM}}$	$1.18$	$1.12$	$1.11$
Computer: CPU Threadripper 1950 ($3.4$ GHz), $128$ GB RAM.

Equations

$$p\left(\mathbf{I}\right)=h\left(\left[\begin{array}[]{*{30}{c}}\cdots,&\mathbf{I}^{\mathrm{H}}\mathbf{A}_{m}\mathbf{I}\,,&\cdots\end{array}\right]\right),$$

$$\tau_{b}\left(p,\varOmega_{\mathcal{E}}\right)=-\Big{(}p\left(\mathbf{I}_{\mathcal{E}}\right)-p\left(\mathbf{I}_{b}\right)\Big{)},$$

$$\left[\mathbf{I}_{{\mathcal{E}}\setminus{\mathcal{B}}_{-}}\right]=\left[{\begin{array}[]{*{40}{c}}\cdots,&\mathbf{y}_{{\mathcal{E}},f}-\displaystyle\frac{Y_{fb}}{Y_{bb}}\mathbf{y}_{{\mathcal{E}},b}\,,&\cdots\end{array}}\right]l_{f}V_{0},$$

$${\mathbf{p}}_{\mathcal{E}\setminus {\mathcal{B}}_{-}}=h\left(\mathcal{D}\left({\mathbf{A}}_{\mathcal{E},m},\left[{\mathbf{I}}_{\mathcal{E}\setminus {\mathcal{B}}_{-}}\right]\right)\right),$$

$$\mathcal{D}\left({\mathbf{A}}_{\mathcal{E}},\left[{\mathbf{I}}_{\mathcal{E}\setminus {\mathcal{B}}_{-}}\right]\right)=\mathrm{diag}\left({\left[{\mathbf{I}}_{\mathcal{E}\setminus {\mathcal{B}}_{-}}\right]}^{\mathrm{H}}{\mathbf{A}}_{\mathcal{E}}\left[{\mathbf{I}}_{\mathcal{E}\setminus {\mathcal{B}}_{-}}\right]\right).$$

$${\mathbf{Y}}_{\mathcal{E}\setminus b}={\mathbf{C}}_{\mathcal{E}\setminus b}^{\mathrm{T}}\left({\mathbf{Y}}_{\mathcal{E}}-\left(\frac{{\mathbf{y}}_{\mathcal{E},b}}{Y_{bb}}\right){\mathbf{y}}_{\mathcal{E},b}^{\mathrm{T}}\right){\mathbf{C}}_{\mathcal{E}\setminus b},$$

$$C_{{\mathcal{E}}\setminus b,nn}=\left\{\begin{array}[]{lll}0&\Leftrightarrow&n=b,\\ 1&\Leftrightarrow&\mathrm{otherwise},\\ \end{array}\right.$$

$$\begin{split}&\left[\mathbf{I}_{{\mathcal{E}}\cup{\mathcal{B}}_{+}}\right]=\\ &\quad\left[\begin{array}[]{*{30}{c}}\cdots,&\mathbf{C}_{{\mathcal{E}}\cup b}^{\mathrm{T}}\left(\left[\begin{array}[]{*{10}{c}}\mathbf{y}_{f}\\ 0\end{array}\right]+\dfrac{x_{fb}}{z_{b}}\left[\begin{array}[]{*{10}{c}}\mathbf{x}_{b}\\ -1\end{array}\right]\right)\,,&\cdots\\ \end{array}\right]l_{f}V_{0},\end{split}$$

$${\mathbf{x}}_b=\mathbf{Y}{\mathbf{z}}_b,z_b=Z_{bb}-{\mathbf{z}}_b^{\mathrm{T}}{\mathbf{x}}_b,$$

$$\mathbf{C}_{{\mathcal{E}}\cup b,mn}=\left\{\begin{array}[]{lll}1&\Leftrightarrow&n=S\left(m\right),\\ 0&\Leftrightarrow&\mathrm{otherwise},\\ \end{array}\right.$$

$$\mathbf{Y}_{{\mathcal{E}}\cup b}=\frac{1}{z_{b}}\mathbf{C}_{{\mathcal{E}}\cup b}^{\mathrm{T}}\left[{\begin{array}[]{cc}z_{b}\mathbf{Y}+\mathbf{x}_{b}\mathbf{x}_{b}^{\mathrm{T}}&\displaystyle-\mathbf{x}_{b}\\ -\mathbf{x}_{b}^{\mathrm{T}}&\displaystyle 1\\ \end{array}}\right]\mathbf{C}_{{\mathcal{E}}\cup b},$$

$$\kappa =\frac{1}{N-1}\Vert {\mathbf{g}}_{\mathrm{init}}\oplus {\mathbf{g}}_{\mathrm{final}}\Vert$$

Full text

I Introduction

Shape synthesis, a technique of constructing a particular body from piece-wise constant materials, is an unsolved problem across many engineering branches. The major obstacle is a $2^{N}$ combinatorial explosion [1] for $N$  degrees-of-freedom (d-o-f) that arises from its binary nature: each unknown is associated with a given material or with vacuum [2]. While the formulation can be relaxed by introducing continuous variables, as in the case of topology optimization [3], the solution is finally rounded with respect to a given threshold [4]. This last step is encumbered with difficulties such as non-uniqueness and instability [5].

Contemporary solutions to shape optimization are mostly parametric sweeps, heuristic algorithms [6], and, recently, machine learning [7]. All these techniques share a common feature: a demand on vast amounts of samples for which a fitness function has to be evaluated. Therefore, large datasets are dealt with during shape optimization.

In this contribution, the binary nature of the optimization problem is kept in its original form accepting NP complexity, and a novel method of topology sensitivity [8] is adopted and extended by the possibility of reconstructing the structure. The resulting local algorithm is based on an investigation of nearest neighbors in a Hamming graph $H(N,2)$ [9]. It utilizes method of moments (MoM) [10] and the Sherman-Morrison-Woodbury identity [11, 12], a popular scheme for evaluating consequences of local geometry perturbations [8, 13, 14, 15, 16].

Exploration of all shapes with a unit Hamming distance is performed with an inversion-free evaluation. The proposed method is fine-grained (with an expected linear speed-up when parallelization is used), well-suited for vectorization, and its implementation has the potential to evaluate millions of mid-size antenna candidates per minute on a laptop. As such, it represents an ideal data mining tool for machine learning algorithms [17] or an apt candidate for a local step in global optimization [18]. Specifically, the training phase of supervised learning [19] can be significantly shortened. Utilizing the linear regression models [19], the proposed technique can provide additional information about first-order perturbations. These claims are supported by a Monte Carlo simulation with a greedy search over the nearest neighbors.

The letter is organized as follows. Topology sensitivity technique is briefly reviewed in Section II and extended by the possibility of reconstructing a previously reduced shape. Section III shows that an iterative evaluation of all nearest neighbours of an actual shape can be employed in a greedy algorithm. Performance of the greedy algorithm based on the nearest neighbours search is statistically studied in Section IV. The letter concludes in Section V.

II Effective Shape Perturbation in MoM

MoM [10] for a fixed discretization and a set of piece-wise basis functions $\left\{\boldsymbol{\psi}_{n}\left(\boldsymbol{r}\right)\right\}$ is considered here as the starting point, see Fig. 1. The presence of basis functions is a subject of binary optimization, see [8] for details. The primary quantities being operated on are the impedance matrix  $\mathbf{Z}\in\mathbb{C}^{N\times N}$ [10], which carries information about the electromagnetic behaviour of an actual shape, and a set of matrix operators $\left\{\mathbf{A}_{m}\right\}$ representing criterion function $p$ in a form

[TABLE]

where $\mathbf{I}$ is a column vector of current expansion coefficients [10], and $h$ is an arbitrary function operating over individual quantities $\mathbf{I}^{\mathrm{H}}\mathbf{A}_{m}\mathbf{I}$, $m\in\left\{1,\dots,M\right\}$. Let us further assume that the $f$th basis function is fed where index $f$ is fixed.

Topology sensitivity is defined as [8]

[TABLE]

expressing a difference between the performance of the actual shape $\varOmega_{\mathcal{E}}$ and of its smallest perturbations realized via individual removals or additions of the $b$th basis function, $b\in\mathcal{B}$. The symbol ${\mathcal{E}}$ denotes a set of all “enabled” d-o-f and ${\mathcal{B}}$ denotes a set of all investigated edges, ${\mathcal{B}}=\left\{{\mathcal{B}}_{-},{\mathcal{B}}_{+}\right\}$, where ${\mathcal{B}}_{-}$ is a set of d-o-f to be removed (one by one) and ${\mathcal{B}}_{+}$ is a set of d-o-f to be added (one by one), see Fig. 2 for a particular example.

II-A Removal of Basis Functions

Shape reduction via the removal of d-o-f with indices $b\in{\mathcal{B}}_{-}$ is possible with a linear asymptotic complexity as

[TABLE]

with an effective evaluation of (1) as

[TABLE]

where

[TABLE]

The symbols introduced in (3)–(5) read: $\mathbf{y}_{{\mathcal{E}},f}$ and $\mathbf{y}_{{\mathcal{E}},b}$ are the $f$th and  $b$th columns of the admittance matrix $\mathbf{Y}_{\mathcal{E}}=\mathbf{Z}_{\mathcal{E}}^{-1}$, respectively, $Y_{fb}$ and $Y_{bb}$ are admittance matrix elements, $l_{f}$ is the length of the fed edge, and $V_{0}=1$ V. Matrix operators $\left\{\mathbf{A}_{{\mathcal{E}},n}\right\}$ have only lines and columns corresponding to the entries in set ${\mathcal{E}}$ and matrices in a form of $\left[\mathbf{I}\right]$ contain columns of current expansion coefficients corresponding to the investigated perturbations.

The edge corresponding to the worst topology sensitivity is selected (here denoted by index $b$) and removed (disabled) by virtue of an admittance matrix update

[TABLE]

where $\mathbf{C}_{{\mathcal{E}}\setminus b}$ is a permutation matrix,

[TABLE]

in which all zero columns are removed. A removal of one basis function thus means reduction of the admittance matrix dimension by one.

II-B Addition of Basis Functions

The basis function removal technique (3) is excellent for investigating topological sensitivity (2), however, a nonexistence of a technique of a basis function addition commonly caused a premature deadlock [8]. Here, the addition of a basis function is introduced by further applying the Sherman-Morrison-Woodbury identity [11] which results in

[TABLE]

where

[TABLE]

and permutation matrix $\mathbf{C}_{{\mathcal{E}}\cup b}$ provides a correct ordering of the basis functions since the basic Sherman-Morrison-Woodbury identity demands that the basis function added must be the last one. Entries of $\mathbf{C}_{{\mathcal{E}}\cup b}$ read

[TABLE]

with $m\in\left\{1,2,\dots,E+1\right\}$ and where $S$ is a set of target indices if a set $\left\{{\mathcal{E}},b\right\}$ is sorted in ascending order.

After deciding which basis function should be added (enabled), an admittance matrix update is performed as follows

[TABLE]

with the auxiliary variables defined in (9).

III Local Shape Perturbation

Thanks to (3)–(6), and (8)–(11) the initial shape $\varOmega_{\mathcal{E}}$ can either be extended or reduced according to its actual topology sensitivity (2) to a given parameter $p$. In order to proceed further, let us represent any properly discretized shape as a binary genus $\mathbf{g}=\left[b_{1},\,\,\cdots\,\,,b_{N}\right]$ with logical values $b_{n}\in\left\{0,1\right\}$, where $b_{n}=1$ if $n\in{\mathcal{E}}$ and $b_{n}=0$ otherwise. The Hamming graph $H(2,N)$ is defined over genes $\mathbf{g}$ in which nearest neighbors can be found and evaluated using (3) and (8), see Fig. 3 for a sketch of the procedure for $N=4$. For an arbitrary starting position, a graph in Fig. 3 can be explored for a locally optimal shape using a greedy algorithm [21] following the steepest descent of (2).

A particular result of the greedy algorithm, based on nearest neighbors, is presented in Fig. 4 for the starting position depicted in Fig. 2 and a minimization of a radiation Q-factor which is a parameter of primary importance for electrically small antenna, [22]. Electrical size is $ka=1/2$, with $k$ being the free-space wavenumber and $a$ being the radius of the sphere fully circumscribing the rectangular region. Actual performance in Q-factor, evaluated according to [23], is normalized throughout the paper with respect to its lower bound being restricted to TM modes only, $Q_{\mathrm{lb}}^{\mathrm{TM}}$, see [24] and references therein. As compared to Fig. 2, the performance of a shape in Fig. 4 was improved from $Q_{\mathrm{init}}/Q_{\mathrm{lb}}^{\mathrm{TM}}=119$ to $Q_{\mathrm{final}}/Q_{\mathrm{lb}}^{\mathrm{TM}}=1.23$ in $84$ steps. The difference between the initial genome $\mathbf{g}_{\mathrm{init}}$ and the final genome $\mathbf{g}_{\mathrm{final}}$ is significant, cf. Figs. 2 and 4, and is quantified with a normalized Hamming distance between the corresponding genes, i.e.,

[TABLE]

for one basis function being excited by a delta gap. The coefficient $\kappa=1$ means the structure was completely changed, while $\kappa=0$ means that the structure was not modified at all. In the final case depicted in Fig. 4, the $\kappa$ coefficient equals $0.47$, which means that approximately one half of the basis functions have been changed when compared with the initial structures.

IV Monte Carlo Analysis

This section provides a detailed study of the algorithm proposed in the previous section. To this point, multiple runs with random starting positions (feeding position being fixed) were performed and statistically evaluated. As in the previous section, the performance of the algorithm is investigated using Q-factor minimization and a rectangular bounding box with a side aspect ratio of 2:1 with electrical size $ka=1/2$. Two discretization schemes, $4\times 8$ grid ($N=180$) and $6\times 12$ grid ($N=414$), are used as can be seen in the insets in Fig. 5 and Fig. 7, respectively. The entire procedure was implemented in MATLAB [26]. The matrix operators were evaluated in AToM [27], and all calculations ran on the computer specified in Table I. The only run-time variables are $\mathbf{g}_{\mathcal{E}}$, $\mathbf{Y}_{\mathcal{E}}$, $\mathbf{Z}$ (for shape reconstruction), and $\left\{\mathbf{A}_{m}\right\}$ matrices. In total, $5\cdot 10^{4}$ trials were performed for both grids and the overall performance is summarized in Table I and in Figs. 5–8.

It can be seen in Fig. 5 that the greedy search, which is a local algorithm, is capable of improving the performance in Q-factor to a mean value of $0.064Q_{\mathrm{init}}$ and $0.028Q_{\mathrm{init}}$ for $4\times 8$ and $6\times 12$ grids, respectively. As expected, the improvements are more pronounced when working with a finer grid (non-hatched bars). As confirmed in Fig. 6, this is only possible if both removal and addition techniques are involved. In particular, the removal alone performs poorly when starting from a random seed. Figure 7 shows the statistics of resemblance between the initial and final samples and reveals that, on average, $\kappa\approx 0.43$ of all d-o-f were modified. The final result, showing the probability density function (PDF) and the cumulative density function (CDF) of a normalized Q-factor, is presented in Fig. 8. Interestingly, the most probable Q-factor is close to its lower bound and is reachable from many starting positions.

V Conclusion

An inversion-free method was introduced for evaluation of the smallest perturbation within the MoM. It makes it possible to preserve the binary nature of the shape synthesis problem, furthermore improving the convergence rate and robustness of the optimization method. The reconstruction of the structure was derived using the Sherman-Morrison-Woodbury identity.

A greedy algorithm used on topology sensitivity was employed to demonstrate the capability to gather millions of evaluated shapes per minute. To this end, the presented method can be utilized as a local step in global optimization schemes. When randomly restarted, it can also serve as a data mining tool or as a building block for machine learning techniques aimed at shape/pattern synthesis.

The letter also stimulates further development. It is inevitable that a study of the dependence of the method on the number of unknowns and the type of mesh grid will be required. Another important question to be discussed concerns the multi-objective formulation and the proposal of a hybrid optimization algorithm based on a combination of a heuristic approach and topology sensitivity.