# Fast Exact Retrieval for Nearest-neighbor Lookup (FERN)

Richard Zhu
Princeton University
<EMAIL_ADDRESS>

###### Abstract

Exact nearest neighbor search is a computationally intensive process, and even
its simpler sibling — vector retrieval — can be computationally complex. This
is exacerbated when retrieving vectors which have high-dimension $d$ relative
to the number of vectors, $N$, in the database. Exact nearest neighbor
retrieval has been generally acknowledged to be a $O(Nd)$ problem with no sub-
linear solutions. Attention has instead shifted towards Approximate Nearest-
Neighbor (ANN) retrieval techniques, many of which have sub-linear or even
logarithmic time complexities. However, if our intuition from binary search
problems (e.g. $d=1$ vector retrieval) carries, there ought to be a way to
retrieve an organized representation of vectors without brute-forcing our way
to a solution. For low dimension (e.g. $d=2$ or $d=3$ cases), kd-trees provide
a $O(d\log N)$ algorithm for retrieval. Unfortunately the algorithm
deteriorates rapidly to a $O(dN)$ solution at high dimensions (e.g. $k=128$),
in practice. We propose a novel algorithm for logarithmic Fast Exact Retrieval
for Nearest-neighbor lookup (FERN), inspired by kd-trees. The algorithm
achieves $O(d\log N)$ look-up with 100% recall on 10 million $d=128$ uniformly
randomly generated vectors.111Code available at
https://github.com/RichardZhu123/ferns

Fast Exact Retrieval for Nearest-neighbor Lookup (FERN)

Richard Zhu Princeton University<EMAIL_ADDRESS>

##  1 Introduction

Vector retrieval is pervasive, underlying search engines, transformers, and
open-book language models. At heart, one of key attributes of computing
systems lie in their ability to retrieve knowledge. Sometimes, when this
knowledge is sufficiently broad — and a powerful enough retrieval architecture
is built — these computing systems may even be so good at retrieving relevant
knowledge that they appear to reason Bubeck et al. (2023).

Given the power of knowledge retrieval for both commercial and academic
pursuits, significant energy has been devoted towards effectively converting
various types of data into vectors, from words Mikolov et al. (2013), images
Radford et al. (2021), and audio Radford et al. (2022) to documents,
sentences, and paragraphs Dai et al. (2015).

In this work, we differentiate between look-up and search. Look-up involves
the retrieval of vectors guaranteed to be contained in the database, while
search involves the retrieval of queries not necessarily contained in the
database. Retrieving queries without exact matches may involve instead
retrieving that vector’s nearest neighbors. The definition of nearest can be
further disambiguated into Euclidean or cosine similarity measures, among
others. Note that under the hood, a Euclidean distance-based nearest neighbor
algorithm can be easily adapted to be cosine similarity-based simply by
dividing each vector in the database by its magnitude during insertion. During
look-up, the query vector is then also divided by magnitude. The resulting
nearest neighbors we obtain are also such by cosine similarity.

The scalable and effective look-up of large numbers of high dimensional
vectors is thus desired. While the vanilla hashmap algorithm provides $O(1)$
time complexity for scalars, extending to $O(d)$ for vectors in d-dimensional
space, this holds only when the cardinality of the hash function range is
large relative to the number of elements, $N$. When the number of elements
becomes large relative to the number of bins, $b$, finding the key within each
bin becomes a linear search problem. Since we expect each bin to have
$\frac{N}{b}$ collisions, the time complexity for look-up is
$O(\frac{Nd}{b})$. While we can also get fast look-up with a heap in low
dimensions, FERN can be thought of as an extension of the heap to high
dimensions while also presenting a novel approach that could lead to sub-
linear vector search.

##  2 Related work

Prior work has taken 3 major directions to attain fast nearest neighbor
search. The approaches involve bucketing, divide and conquer, or graph-based
approaches. These techniques are specifically tuned to work well for vectors
in high dimensional space, which should be unsurprising since each of these
techniques is really just an extension of familiar 1-d concepts: hashmaps,
binary search, and breadth first search. I believe there are key learnings
that can be taken from both exact and approximate retrieval and search
settings, even though we are interested in the exact variant. Some new
techniques such as Certified Cosine certificates Francis-Landau and Van Durme
(2019) offer structure-agnostic tweaks to speed up performance.

###  2.1 Bucketing

Locality Sensitive Hashing Indyk and Motwani (1998) and k-means Lloyd (1982);
MacQueen (1967), including more recent variants like k-means++ Arthur and
Vassilvitskii (2007), are two techniques that use bucketing to group database
vectors with their nearby peers. Since it is hard to deal with queries that
fall on the boundaries of adjacent clusters, these algorithms are approximate
search algorithms. Both of these algorithms in practice take linear time for
sufficiently high dimensions and large numbers of elements.

#### Divide and conquer

kd-trees provide strong algorithms for pruning, with new variants attempting
to decrease constant factors in the time complexity. Zhang et al. (2012) A
recent work Ram and Sinha (2019) attempts to create partitions based on random
rotations of the dataset achieve the same search accuracy guarantees as RPTree
Dasgupta and Sinha (2013) but with $O(d\log d+\log n)$ time complexity for
approximate search.

#### Graph-based approaches

Recent graph-based approaches to search, such as Navigable Small World (NSW)
graphs Malkov et al. (2014) and a later variant involving layers of NSWs,
provide increasing granularity. NSWs are graphical representations of
databases where each pair of nodes is connected by a small number of hops.
Further optimizations have been presented Fu et al. (2019); Jayaram Subramanya
et al. (2019).

##  3 FERN

#### Goal

We aim to build an algorithm that satisfies two primary goals: we must be able
to perform quick look-up on vectors guaranteed to be contained in the database
and we must be able to quickly insert vectors. When designing the algorithm,
we initially assume our database contains $N$ vectors, each spanning $d$
dimensional space. The $i$-th vector $v_{i}$ is defined as follows

$\displaystyle\mathbf{v_{i}}=[v_{i,1},v_{i,2},\ldots,v_{i,d}]^{\top}$
$\displaystyle\text{where }v_{i}\sim\mathcal{R}(-1,1)$

This process effectively generates vectors lying within a $d$ dimensional
"ball"-like shape of radius 1 in each direction. The directions in which the
vectors point are also evenly distributed direction-wise. We discover later
however, that the algorithm we design with this simplifying assumption
maintains logarithmic time lookup for $\mathbf{v_{i}}$ of any arbitrary
direction and length.

In terms of time complexity, we define quick as anything taking logarithmic
time — this means that lookup in a database of ten billion vectors should only
take seven times longer than lookup in a database of a thousand vectors. That
is remarkable, because a naive linear search would take ten million times
longer, a nearly intractable time difference at scale. Since a logarithmic
lookup time and linear space complexity is state of the art (SoTA), we believe
a key contribution of our work is an alternative data structure and algorithm
that achieves SoTA while simultaneously providing the capacity to be extended
to logarithmic-time nearest neighbor search, given its unique approach to
dividing the vectors by hyperplanes defined based on the vectors in the
database rather than measuring along a specific direction like the traditional
kd-trees process. Each node is an object that stores a vector, pointers to the
left and right children, and a pointer to that node’s parent node. This
results in a binary tree with undirected edges. While we ultimately implement
retrieval using a queue structure, this bidirectional edge only adds marginal
complexity to the algorithm and underlying data structure while enabling a
backtracking-based traversal method. The queue-based method emulates a level-
order traversal of candidate nodes while a stack-based backtracking-based
traversal method (that fully explores a specific path before backtracking;
exploring each sibling node that could not be pruned without potentially
missing the nearest neighbor) emulates a depth-first search.

#### Methodology

We design a novel algorithm that is a variant of kd-trees, but has the
capacity to perform better at higher dimensions. Broadly, the structure is a
binary tree. Each node that has both left and right children defines a
hyperplane using the vectors of its left and right children as support
vectors.

The tree is constructed so that all vectors in each child’s subtree are on the
same side of the hyperplane as that child. This allows us to perform vector
look-up in logarithmic time, provided that vectors are added to the database
in a sufficiently random manner. It is feasible, however, for an adversarial
insertion process to result in a heavily imbalanced tree and consequently
worst-case linear look-up time. While we don’t observe this as an issue in
practice, we can resolve the issue by implementing a slightly more complicated
variant of FERN - using a variant of the Red-Black Tree technique Guibas and
Sedgewick (1978) to guarantee balanced trees, logarithmic depth, and thus
logarithmic retrieval time complexity.

There are two key components to our algorithm: one method for insertion
(Algorithm 1) and another for lookup (Algorithm 2).

The insertion algorithm (Algorithm 1) is fairly concise. When inserting a
vector into the tree, it is placed at the root if the tree has not been
initialized yet. Otherwise, if the current node is missing a left or right
child, we insert the vector as a child node. If the node is a leaf node — that
is, missing both left and right children — then the left child is always
inserted first, before the right child.

Algorithm 1 FERN Insertion

1: Function Insert(vector)

2: if root not initialized then

3: $\text{root}\leftarrow\text{VectorNode(null)}$

4: end if

5: $\text{node}\leftarrow\text{root}$

6: while True do

7: if no left child then

8: set left node to vector

9: break

10: else if no right child then

11: set right node to vector

12: break

13: else if vector closer to left child then

14: $\text{node}\leftarrow\text{node.left}$

15: else

16: $\text{node}\leftarrow\text{node.right}$

17: end if

18: end while

During insertion, if a node has both left and right children, then we set the
current node instead to the child node that is closest to the vector we are
inserting. That is, if we form a hyperplane from the set of points equidistant
to both left and right children, then we set the current node to the left
child if the vector to be inserted lies on the same side of the hyperplane as
the left child, otherwise we set the current node to the right child.

The result of this insertion algorithm is that — for a balanced binary tree —
we get a maximum tree depth of $O(\log N)$ where $N$ is the number of elements
in the database. Insertion time per vector is thus $O(\log N)$ since we only
visit one node per depth level.

When looking up vectors from the data structure, we demonstrate a method
(Algorithm 2) that has a per-vector retrieval time proportional to the maximum
depth of the tree, since we only look at one node per depth level. However,
when extended to search settings where the query is not known to be contained
in the database, retrieval time becomes proportional to the number of elements
in the tree. We can no longer automatically prune any queries that lie close
to the hyperplane boundary since there is a possibility that the nearest
neighbor and query may lie on different sides of the hyperplane.

Ostensibly, we would expect a non-negligible proportion of vectors to be
sufficiently far from the hyperplane to be pruned. However, we notice that in
practice, as the dimensionality of the vectors increase, so too does the
proportion of vectors lying close to the boundary. This makes sense
intuitively since we are trying to project increasingly higher dimensions of
vectors onto a 1-d line (the line normal to the hyperplane and passing through
both support vectors). During retrieval, we effectively perform the depth-
first or level-order search, as described previously. For a balanced tree with
strong boundaries (that is, most queries lie far away from the hyperplane),
per-vector time should be logarithmic with respect to the number of elements
already present in the database (another word for our proposed data structure)
at insertion-time. However, it becomes linear otherwise. We first define mip
and mip_vec, the distance to the nearest neighbor found thus far and the
vector representing the nearest vector retrieved thus far. We then create a
queue and insert the root node.

Algorithm 2 FERN Lookup

0: query

0: mip_vec

1: $\text{mip},\text{mip\\_vec}\leftarrow\infty,\text{None}$

2: $\text{curr}\leftarrow\text{None}$

3: $\text{queue}\leftarrow[\text{self.root}]$

4: while queue not empty do

5: $\text{curr}\leftarrow\text{oldest element in queue}$

6: update mip, mip_vec if curr closer to query

7: if curr has both left and right children then

8: if query is closer to left child then

9: queue.append(curr.left)

10: if query close to boundary then

11: queue.append(curr.right)

12: end if

13: else

14: queue.append(curr.right)

15: if query close to boundary then

16: queue.append(curr.left)

17: end if

18: end if

19: else if curr has left child only then

20: queue.append(curr.left)

21: else if curr has right child only then

22: queue.append(curr.right)

23: end if

24: end while

25: return mip_vec

We then continuously pop a node from the head of the queue until the queue is
empty. Each time we pop a node, we check whether its vector is closer to the
query vector than the current best candidate for nearest neighbor, mip_vec,
which is a Euclidean distance mip away from the query. For lookup, we are
looking for an exact match, so we are seeking an mip of 0. If the node has a
left child only or a right child only (the latter should never happen, but we
have it as a redundancy against exceptions) then we add that node to the
queue. Otherwise, if both children exist then we add to the queue the node
that shares the same side of the hyperplane as the query. For lookup cases, we
consider any query to be sufficiently far from the boundary that only one
child node needs to be added to the queue per node. After all, whether a query
is close to the boundary is somewhat arbitrary and the exact function
definition depends on whether we are performing lookup or search.

For mapping applications, we can add an additional variable, data (a byte
array), to the Node class.

In Algorithms 1 and 2, “closeness" is quantified by Euler distance, and the
“boundary" is the midpoint between left and right child nodes.

##  4 Experimental results

During experiments, we typically utilize $d=128$ with the same uniform
distribution previously assumed. While this may not be characteristic of all
data distributions, we note that our architecture is actually agnostic to the
distribution of the vectors being inserted. What matters (in terms of
potential effects on performance) is the order in which vectors are inserted
based on their relative positions.

To run our experiments, we use the Intel Xeon Platinum 8380 CPU (2.30 GHz),
the same processor used for running the popular ann-benchmark Bernhardsson
(2024). For values of $N$ we use $10^{4}$, $5*10^{4}$, $10^{5}$, $5*10^{5}$,
and $10^{6}$, $5*10^{6}$, $10^{7}$. The last setting, equivalent to look-up on
10 million vectors, has comparable values of $N$ and $d$ to many of the
Euclidean distance based benchmarks in ann-benchmark. We notice a nearly
perfect logarithmic time complexity, and at $N=10^{7}$ we run approximately
3000 retrievals/second without additional optimizations.

##  5 Discussion

During our experiments, we noticed that in the search modality, system
dynamics can change drastically based on the dimensionality of the vectors. We
experiment with different ways of defining the hyperplanes and various
algorithms that would balance or repair the tree to try to guarantee
logarithmic retrieval for large $d$ and $N$. Unfortunately, while these
algorithms almost universally gave logarithmic time complexity for 100%
recall, the performance broke down drastically beyond $d=2$ or $d=3$. In
particular, we note the importance of having well-defined hyperplane
boundaries.

#### Boundary sharpness

We want to be able to maximize pruning since we achieve $O(\log_{2}N)$ time
complexity when we prune 50% of the nodes in the database each time we measure
the distance between a node and the query). Indeed, note that since a
$k$-means based nearest neighbors search allows us to prune up to $N/k$ nodes
per comparison, we might wonder why $2$-means search doesn’t achieve $O(\log
N)$ complexity. It’s because $k$-means has an $O(N/k+k)$ time complexity, the
$N/k$ term means that we would need to do a linear number of searches
regardless of the number of clusters.

We observe empirically that the proportion of nodes in the database that are
visited increases sharply when there are more nodes that are closer to the
hyperplane than the support vectors that define the plane. That is to say,
during the insertion process, there will be vectors that lie between a support
vector and the hyperplane. Now when we’re retrieving, the query may lie on one
side of the hyperplane but its nearest neighbor may be one of these "in-
between" nodes on the other side of the hyperplane. This means we now must be
much more prudent when pruning which decreases the proportion of vectors that
are pruned and thus increases the time complexity.

However in the look-up modality, we achieve logarithmic time complexity on
both a vector database of dimension $d=128$ and size $10^{7}$ that are
randomly and uniformly generated (Figure 1) and three Euclidean benchmarks
(Table 1) from ann-benchmark (Figure 2), using the larger train splits
($N=60,000$ to $N=1,000,000$). We observe logarithmic time complexity over a
diverse dimensions and vector distributions.

Dataset | Dim. | Train Size | Test Size
---|---|---|---
Fashion-MNIST | 784 | 60,000 | 10,000
MNIST | 784 | 60,000 | 10,000
SIFT | 128 | 1,000,000 | 10,000
Table 1: Properties of ANN Benchmark datasets used for evaluation Figure 1:
FERN lookup with vectors where $d=128$ and look-up time is averaged over 1000
vectors randomly sampled from the database Figure 2: FERN lookup using the
train portion (60k-100k vectors) of popular Euclidean-distance-based vector
retrieval benchmarks and look-up time is averaged over 1000 vectors randomly
sampled from the database. We evaluate 4 decades on each dataset, which is why
SIFT-128-Euclidean evaluation starts with vector databases of size $10^{3}$
rather than $600$

##  6 Conclusion

We are able to achieve our goal of creating a novel vector database structure
that achieves state of the art look-up time complexity that is logarithmic in
the number of vectors. The algorithm presented here, FERN, maintains 100%
recall while performing lookup on vectors in high-dimensional space (e.g.
$d=128$ to $d=784$) with $N$ varying from $10^{2}$ to $10^{7}$, and presents a
potential path towards a data structure and algorithm that will allow for the
first sub-linear exact nearest neighbor retrieval process.

We believe that the exact process for attempting to perform binary search on a
vector database requires carefully defined hyperplanes, which presents an area
for further work. We find the "fixing" step of the Red-Black tree algorithm to
be particularly inspirational as a direction of future work. Evaluation on
additional datasets and further investigations of recall-retrieval-time trade-
offs could help in the pursuit of sub-linear search.

We further also believe that an alternative for hyperplanes is to use a graph
based approach, similar to the approach taken in many recent works Malkov et
al. (2014); Fu et al. (2019); Jayaram Subramanya et al. (2019), since this
could allow us to more easily divide the database in a well-defined and easy
to update way. Overall we are excited by the potential and hope to further
develop this algorithm in pursuit of sub-linear exact search.

## Acknowledgements

The author thanks Edo Liberty, Huacheng Yu, Runzhe Yang, Nataly Brukhim,
Douglas Downey, and Matthijs Douze for their invaluable advice and feedback
during the writing of this paper.

## References

  * Arthur and Vassilvitskii (2007) David Arthur and Sergei Vassilvitskii. 2007.  k-means++: the advantages of careful seeding.  In _Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms_ , SODA ’07, pages 1027–1035, USA. Society for Industrial and Applied Mathematics.
  * Bernhardsson (2024) Erik Bernhardsson. 2024.  erikbern/ann-benchmarks.  Original-date: 2015-05-28T13:21:43Z.
  * Bubeck et al. (2023) Sébastien Bubeck, Varun Chandrasekaran, Ronen Eldan, Johannes Gehrke, Eric Horvitz, Ece Kamar, Peter Lee, Yin Tat Lee, Yuanzhi Li, Scott Lundberg, Harsha Nori, Hamid Palangi, Marco Tulio Ribeiro, and Yi Zhang. 2023.  Sparks of artificial general intelligence: Early experiments with gpt-4.
  * Dai et al. (2015) Andrew M. Dai, Christopher Olah, and Quoc V. Le. 2015.  Document Embedding with Paragraph Vectors.  ArXiv:1507.07998 [cs].
  * Dasgupta and Sinha (2013) S. Dasgupta and Kaushik Sinha. 2013.  Randomized partition trees for exact nearest neighbor search.  _ArXiv_.
  * Francis-Landau and Van Durme (2019) Matthew Francis-Landau and Benjamin Van Durme. 2019.  Exact and/or Fast Nearest Neighbors.  ArXiv:1910.02478 [cs].
  * Fu et al. (2019) Cong Fu, Chao Xiang, Changxu Wang, and Deng Cai. 2019.  Fast approximate nearest neighbor search with the navigating spreading-out graph.  _Proceedings of the VLDB Endowment_ , 12(5):461–474.
  * Guibas and Sedgewick (1978) Leo J. Guibas and Robert Sedgewick. 1978.  A dichromatic framework for balanced trees.  In _19th Annual Symposium on Foundations of Computer Science (sfcs 1978)_ , pages 8–21.  ISSN: 0272-5428.
  * Indyk and Motwani (1998) Piotr Indyk and Rajeev Motwani. 1998.  Approximate nearest neighbors: towards removing the curse of dimensionality.  In _Proceedings of the thirtieth annual ACM symposium on Theory of computing - STOC ’98_ , pages 604–613, Dallas, Texas, United States. ACM Press.
  * Jayaram Subramanya et al. (2019) Suhas Jayaram Subramanya, Fnu Devvrit, Harsha Vardhan Simhadri, Ravishankar Krishnawamy, and Rohan Kadekodi. 2019.  DiskANN: Fast Accurate Billion-point Nearest Neighbor Search on a Single Node.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc.
  * Lloyd (1982) S. Lloyd. 1982.  Least squares quantization in PCM.  _IEEE Transactions on Information Theory_ , 28(2):129–137.  Conference Name: IEEE Transactions on Information Theory.
  * MacQueen (1967) J. MacQueen. 1967.  Some methods for classification and analysis of multivariate observations.
  * Malkov et al. (2014) Yury Malkov, Alexander Ponomarenko, Andrey Logvinov, and Vladimir Krylov. 2014.  Approximate nearest neighbor algorithm based on navigable small world graphs.  _Information Systems_ , 45:61–68.
  * Mikolov et al. (2013) Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean. 2013.  Efficient Estimation of Word Representations in Vector Space.  ArXiv:1301.3781 [cs].
  * Radford et al. (2021) Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, Gretchen Krueger, and Ilya Sutskever. 2021.  Learning Transferable Visual Models From Natural Language Supervision.  ArXiv:2103.00020 [cs].
  * Radford et al. (2022) Alec Radford, Jong Wook Kim, Tao Xu, Greg Brockman, Christine McLeavey, and Ilya Sutskever. 2022.  Robust Speech Recognition via Large-Scale Weak Supervision.  ArXiv:2212.04356 [cs, eess].
  * Ram and Sinha (2019) Parikshit Ram and Kaushik Sinha. 2019.  Revisiting kd-tree for Nearest Neighbor Search.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 1378–1388, Anchorage AK USA. ACM.
  * Zhang et al. (2012) Yuezhuo Zhang, Li Zha, Jia Liu, and Zhiwei Xu. 2012.  A Large-Scale Online Search System of High-Dimensional Vectors Based on Key-Value Store.  In _2012 Eighth International Conference on Semantics, Knowledge and Grids_ , pages 233–236.


11institutetext: National Solar Observatory, 3665 Discovery Drive, 3rd Floor,
Boulder, CO 80303 USA
11email<EMAIL_ADDRESS>22institutetext: Stanford University, Stanford, CA,
USA 33institutetext: ReSoLVE Centre of Excellence, Space Climate research
unit, University of Oulu, POB 3000, FIN-90014, Oulu, Finland 44institutetext:
NorthWest Research Associates, 3380 Mitchell Lane, Boulder, CO 80301 USA
55institutetext: Institute for Space-Earth Environmental Research, Nagoya
University, Furo-cho Chikusa-ku Nagoya, Aichi 464-8601 Japan

# On a limitation of Zeeman polarimetry and imperfect instrumentation in
representing solar magnetic fields with weaker polarization signal

Alexei A. Pevtsov 11 Yang Liu 22 Ilpo Virtanen 33 Luca Bertello 11 Kalevi
Mursula 33 K.D. Leka 44 5 5 Anna L.H. Hughes 11

###### Abstract

Full disk vector magnetic fields are used widely for developing better
understanding of large-scale structure, morphology, and patterns of the solar
magnetic field. The data are also important for modeling various solar
phenomena. However, observations of vector magnetic fields have one important
limitation that may affect the determination of the true magnetic field
orientation. This limitation stems from our ability to interpret the differing
character of the Zeeman polarization signals which arise from the photospheric
line-of-sight vs. the transverse components of the solar vector magnetic
field, and is likely exacerbated by unresolved structure (non-unity fill
fraction) as well as the disambiguation of the 180∘ degeneracy in the
transverse-field azimuth. Here we provide a description of this phenomenon,
and discuss issues, which require additional investigation.

###### keywords:

Sun: magnetic fields – Techniques: polarimetric

##  1 Introduction

Observations of magnetic fields provide key information for developing our
understanding of the Sun’s short-term (space weather) and long-term (space
climate) activity and in predicting these effects on Earth. Synoptic full disk
longitudinal magnetograms have existed since the late 1960s, and these data
continue to serve as the primary input for space weather and space climate
research and operational forecasts. By their nature, longitudinal magnetograms
do not contain sufficient information to derive the true orientations of the
magnetic-field vectors, and thus, require additional assumptions for physical
interpretation. For example, “radial field” synoptic maps, which are widely
used in space weather forecasting are created under the assumption that the
true field is radial. Observations of vector Stokes polarimetery in principle
have the information necessary to fully reconstruct photospheric vector-
magnetic-field maps, and efforts are underway to employ such data in
operational space weather forecasts.

The earliest observations of vector magnetic fields in solar active regions
were conducted at the Crimean Astrophysical Observatory in the early 1960s
(Stepanov and Severny, 1962; Severny, 1965). By the early 1980s, a number of
vector magnetographs were developed around the world, with the most prolific
instruments operating in Czechoslovakia, East Germany (Pflug and Grigoryev,
1986), Japan (NAOJ, Ichimoto et al., 1993), the Soviet Union (Crimean, Pulkovo
and Sayan observatories), and the USA (NASA’s Marshall Space Flight
Center/MSFC, Mees Solar Observatory of University of Hawaii, High Altitude
Observatory/HAO) (for review, see individual articles in Hagyard, 1985;
Cacciani et al., 1990). All of these instruments had a limited field of view,
typically about the size of an average active region. Full disk vector
magnetograms have been routinely observed since late 2003 by the Vector Stokes
Magnetograph (VSM) on the Synoptic Optical Long-term Investigation of the Sun
(SOLIS) platform (Keller et al., 2003). Beginning in 2010, full disk vector
magnetograms are available from the Helioseismic and Magnetic Imager (HMI,
Scherrer et al., 2012) on board the Solar Dynamics Observatory (SDO).

On 23–26 January 2017, a working meeting on the “Use of Vector Synoptic Maps
for Modeling” was held in Oulu, Finland with two follow up meetings at the
National Solar Observatory, Boulder, Colorado, USA (7–10 November 2017), and
at the Max Planck Institute for Solar System Research, Göttingen, Germany
(18-21 September, 2018). At the first meeting, a direct comparison of vector
magnetic field observations from SDO/HMI and SOLIS/VSM revealed disagreements
in the orientation of meridional (Bθ, North-South) and/or zonal (Bφ, East-
West) components of the vector field in plage areas. One example of this
disagreement is shown in Figure 1, which shows the three components of the
magnetic field in VSM and HMI synoptic maps between heliographic latitudes
-40∘ and 0∘ and Carrington longitudes 60∘ and 130∘.

Figure 1: Example of disagreement in orientation of Bθ and Bφ Magnetic field
vector components of an active region in the southern hemisphere, which passed
the central meridian on 19 November 2015. Upper panels show radial (Br),
meridional (Bθ) and zonal (Bφ) field from SDO/HMI and lower panels from
SOLIS/VSM. Meridional (middle) and zonal (right) components depict opposite
orientation. Red/blue represent positive/negative polarities scaled between
$\pm$ 50 G (see color bars on the right side of figure).

This decaying active region (plage) was at the central meridian on 19 November
2015. A synoptic map is used as an example because the averaging used in
producing a synoptic map results in lower background noise, but very similar
patterns are also seen in less-averaged full-disk data. Recent comparisons
also indicate that the average large-scale zonal field Bφ observed by
SOLIS/VSM and derived from line-of-sight observations tends to disagree
outside active regions (Virtanen et al., 2019). The investigation of this
disagreement uncovered what we think may be an important limitation to Zeeman
vector polarimetry. Here we present the results of our investigation into the
discrepancies uncovered by that collaborative effort. In Section 2, we
introduce observations from two vector magnetographs and compare the vector
field derived from these instruments and their standard data reduction
algorithms. Sections 3-4 introduce our explanation of the observed
disagreement, and in Section 5 we discuss the results of our findings.

##  2 Vector magnetograms from different instruments and a comparative
analysis

Here, we employ full disk vector magnetograms from two instruments: the Vector
Stokes Magnetograph (VSM) on the Synoptic Optical Long-term Investigations of
the Sun (SOLIS) platform (Keller et al., 2003; Balasubramaniam and Pevtsov,
2011), and the Helioseismic and Magnetic Imager (HMI, Scherrer et al., 2012;
Hoeksema et al., 2014) on board the Solar Dynamics Observatory (SDO, Pesnell
et al., 2012).

The HMI instrument is a filtergraph covering the full solar disk with
$4096\times 4096$ pixels. The spatial resolution is about 1” with a 0.5” pixel
size. The width of the filter profiles is 7.6 pm. The spectral line used is Fe
i 617.3 nm, which forms in the photosphere (Norton et al., 2006). The Stokes
parameters $(I,Q,U,V)$ are computed from those measurements (Couvidat et al.,
2016), and are further inverted to retrieve the vector magnetic field using a
Milne-Eddington (ME) based inversion algorithm, the Very Fast Inversion of the
Stokes Vector (VFISV Borrero et al., 2011; Centeno et al., 2014). To suppress
p-mode oscillations and to increase the signal-to-noise ratio, registered
filtergrams are averaged over a certain time before computing the Stokes
vector. Currently a weighted average is computed every 720 seconds using data
obtained over 1350 seconds by default; other averaging windows are available.

Inversion of the vector field has an unavoidable 180∘ ambiguity in the
azimuthal field direction. Assumptions about the field must be made to resolve
the ambiguity. For all pixels in active regions, as well as for strong-field
pixels (where the S/N $>$3 in the transverse signal plus a 50 G buffer) in
quiet Sun regions, the azimuth is determined using a minimum energy algorithm
(Metcalf, 1994; Metcalf et al., 2006; Leka et al., 2009; Hoeksema et al.,
2014). The minimum-energy-method computation is time consuming for pixels
where the signal is dominated by noise, so for weaker polarization regions,
the $180^{\circ}$ ambiguity is solved using three quicker methods: a
randomizing method (the option to add 180∘ is determined randomly), an acute-
angle comparison to a potential field, and a method that provides the most
radially-directed solution. More details can be found in Hoeksema et al.
(2014). In this study, we used the random disambiguation for pixels with
weaker linear polarization.

The VSM is a spectrograph-based instrument, which observes full line profiles
of the Fe i 630.15 and 630.25 nm spectral lines, with a spectral sampling of
2.4 pm and pixel size of 1.0 $\times$ 1.0 (1.14 $\times$ 1.14 before January
2010) arcseconds over a 2048 $\times$ 2048 pixels field of view. To construct
a full-disk magnetogram, the image of the Sun is scanned in the direction
perpendicular to the spectrograph slit. At each scanning step, the spectra for
each pixel along the slit are recorded simultaneously. Each scanning step
takes 0.6 s, and a full disk magnetogram can be completed in about 20 minutes.
The spectrograph slit is curved, which introduces geometric distortions to the
image of the Sun. These distortions are corrected by shifting the position of
each pixel in the final image to the closest integer position of the true
pixel location in a round-Sun image. The maximum uncertainty in the position
of a pixel does not exceed half-a-pixel, which is significantly smaller than
typical atmospheric seeing for this groundbased instrument. The above
correction procedure avoids ill-posed interpolation of full disk magnetograms,
and it preserves the mapping of spectral information for each image pixel.

Similar to HMI, the observed profiles of Stokes Q, U, V, and I parameters are
inverted using the VFISV code under the assumption of a standard Milne-
Eddington stellar atmosphere. However, unlike HMI, VSM inversion includes the
magnetic-field filling factor ($\alpha$) as an additional fit parameter, which
represents the fraction of each instrument pixel filled by magnetized plasma.
For additional details about SOLIS/VSM inversion methods and pipeline, see
Harker (2017). The 180∘ azimuthal ambiguity in the transverse field is
resolved using the Very Fast Disambiguation Method (VFDM, Rudenko and
Anfinogentov, 2014). The VFDM has an accuracy almost as good as that of the
minimum energy method (used for HMI disambiguation), but is much faster. For a
synoptic instrument such as VSM and HMI, the disambiguation is done
automatically as part of the pipeline data reduction. The pipeline reductions
are optimized for “a good answer most of the time, in time for the next
dataset”, and in some cases may not return the best possible solution.

VSM data are from data processing level PROVER0 = 15.0511, which uses only one
(Fe i 630.15 nm) spectral line for the inversion. In the following discussion
we adopt the righthanded coordinate system (Br,Bθ,Bφ), which has been used in
previous publications (e.g., Virtanen et al., 2019). Here the radial component
(Br) is positive when pointing away from the Sun, the meridional component
(Bθ) is positive southward, and the zonal component (Bφ) is positive westward.

##  3 Effects of noise in the transverse fields on the derived vector-field
orientation

Let us now consider how properties of noise could affect the derived
interpretation of the vector magnetic field. In all modern instruments, the
derivation of vector magnetic fields is based on observations of full Stokes
profiles in a selected spectral line sensitive to the effects of the magnetic
field at the location of the formation of this line. Stokes I represents the
total intensity of light. Stokes V is circular polarization (counter-/clock-
wise), and Stokes Q and U represent two linear polarizations. The observed
Stokes profiles are fitted by a model/synthetic line profiles in a process
called “inversion”, and the properties of the magnetic field (and other
parameters, such as Doppler velocity, temperature, magnetic filling factor
etc) are determined based on properties of the fitted line profiles.

The observed profiles of all Stokes parameters are affected by noise and other
instrumental limitations. The photometric noise in Stokes Q and U is similar
to that of Stokes V. However, the longitudinal (B||) field is related linearly
to circular polarization, while the relation between transverse field and
linear polarization is quadratic (e.g., see Equation 18 in Stenflo 2013 and
Equations 12c and 12d in Jefferies and Mickey 1991), which results in a lower
signal/noise in the latter for the same underlying field strength. For
example, in the weak-field approximation (Ronan et al., 1987; Jefferies et
al., 1989; Jefferies and Mickey, 1991),

$B_{||}\approx\frac{C_{||}}{\alpha}~{}{\rm V}$ (1)
$B_{\perp}\approx\frac{C_{\perp}}{\sqrt{\alpha}}\sqrt[4]{{\rm Q}^{2}+{\rm
U}^{2}}.$ (2)

Based on radiative transfer modeling of the Fe i 630.15 nm and assuming a
filling factor of unity, $\alpha=1$ (the entire pixel is filled with a field),
Stenflo (2013) estimated the coefficients of proportionality as $C_{||}\approx
29.4$ and $C_{\perp}\approx 184$. (The coefficients are model dependent, but
$C_{\perp}>>C_{||}$ independent of a model).

As a result, the noise measured in horizontal fields is typically larger than
the amplitude of noise in the longitudinal field by a factor of 10-25.
Moreover, unlike noise in B||, which is distributed around zero,
(B${}_{\perp}\pm$ noise) is always positive (e.g. see Appendix A in Leka and
Skumanich, 1999). For some specific magnetic configurations, this dichotomy
may systematically distort the true inclination of the vector magnetic field.

Figure 2: Effect of noise in the transverse field on the derived orientation
of the vector field. (a) Black arrows represent the true vector field on the
Sun (radial in this example). Red arrows represent the line-of-sight component
and thick blue arrows are transverse components of the true vector field.
Adding positive noise to the transverse field (thin blue arrows) makes the
observed field (dashed arrow) systematically inclined in the direction away
from the central meridian. Because of the 180∘ azimuthal ambiguity in the
transverse field, the same transverse field may satisfy an alternative
orientation (dashed-dotted arrow). The resulting orientation of the field
vector is shown by a dashed arrow for both locations. Letters E, W, N mark
approximate positions of solar East limb, West limb, and North as seen from
Earth. (b) Difference between the true radial direction and vector orientation
in the presence of noise, when azimuth disambiguation selects the most radial
of two possible solutions (open circles). Solid line is a 6th degree
polynomial fit.

Figure 2a shows a theoretical example, where the true magnetic field is
expected to be radial (black arrow). Observed at disc center, a radial
magnetic field is then determined by the observed line-of-sight component plus
symmetric random noise in the azimuth direction, which is negligible on
average. However, outside the central meridian as depicted in Figure 2a, a
purely radial field contributes to both the line-of-sight component (red
arrow) and the transverse component (blue arrow).

In the presence of noise, the orientation of the magnetic vector will be
determined by the observed (B${}_{||}\pm\delta B_{||}$) and
(B${}_{\perp}\pm\delta B_{\perp}$) components, where $\delta B$ represents the
amplitude of noise. Note that due to a quadratic contribution of Stokes Q and
U to B⟂ (see, Equation 2), (B${}_{\perp}\pm\delta B_{\perp}$) is always
positive (by definition, transverse field is always positive, with or without
noise).

For pixels situated East of the central meridian (upper part of Figure 2a),
the projection of the same (radial) vector would create a systematic non-zero
transverse component, which in the presence of noise will only increase (due
to B${}_{\perp}\pm\delta B_{\perp}>$ 0). Due to the 180-degree ambiguity in
the azimuth of the transverse field, the vector field would have two possible
orientations, and if the more radial (black solid arrow) option is chosen,
then the selected orientation of the magnetic field will have a systematic
tilt in the direction away from the solar disk center. Figure 2b shows a
difference between the true radial field and one with $\delta B_{\perp}$ about
20 G added to the true transverse field. For this modeling exercise, we
assumed that the true field strength of the vector field was 200 G, and
$\delta B_{||}$ = 0 (no error in the longitudinal field), and $\delta
B_{\perp}$ was randomly generated within the range of about 0–100 G, with a
mean about 50 G and a standard deviation of 20 G. For pixels located closer to
the disc center, we see an increase in scatter in the vector inclination
relative to the true radial direction. Most importantly, there is a systematic
inclination of the vector field in the direction away from the central
meridian. For a field strength of 200 G and 20 G noise in transverse field,
the inclination error could be up to 20 degrees (Figure 2b).

Figure 3: (Left) Changes in the pattern of the zonal (East-West) vector-field
component of a small bipolar region over its disk passage. Four panels show
the so-called Near Real Time (NRT) synoptic maps. White/black halftones
correspond to magnetic field directed towards the West/East. Each NRT map
covers 360 degrees in longitude (horizontal direction) and approximately $\pm$
30 degrees in latitude (vertical direction). The most recent data are added
onto the left side of synoptic maps. The dates of the most recent observations
added to synoptic map are shown in the upper-right corner of each panel. For
visual clarity, the synoptic charts are shifted to have the active regions
aligned in the vertical direction. Small red arrows plotted in the left part
of each panel correspond to the approximate location of the central meridian
for the day of observations. Panels on the right show a schematic inclination
pattern of the magnetic field vectors in a vertical East-West oriented plane
(as if we were looking at an active region from the side).

Figure 3 demonstrates this effect in the disk-passage of a small bipolar
region with negative leading and positive following polarity fluxes. It is
clear that when this small active region is located to the East of central
meridian, the horizontal magnetic field connecting the two polarities is
directed westward in the leading (negative) polarity and eastward in the
following (positive) polarity flux element. Only when the region is near the
central meridian does it become clear that the magnetic field in both polarity
fluxes is close to radial, with a very slight inclination in the direction
away from each polarity. When the region is located West of central meridian,
the pattern of the zonal field reverses as compared with its location East of
central meridian. This behaviour is in perfect qualitative agreement with the
explanation given in Figure 2.

Let’s now consider whether lowering the noise level will mitigate this
systematic effect. For this test, we selected HMI observations taken on 10
February 2015. That day had a good representation of various solar features
(plage, sunspots) situated at different distances from solar disk center. For
the test, we employed filtergrams that were processed using a much longer
integration than normal (pipeline) magnetograms (5760 s vs. 720 s, courtesy
Dr. X. Sun). Integrated (average) filtergrams then were inverted using
standard VFISV code, and disambiguated using the minimum energy disambiguation
method with default pipeline settings. The averaged magnetogram is centered at
19:12:00 UT, and is shown in Figure 4 (left). Despite significant time
averaging, we do not see any obvious smearing of solar features. The S/N is
improved when compared to the standard 720s magnetograms, especially in areas
with a weak polarization signal.

To test the effect of lower noise on the inferred direction of the magnetic
field vector in large-scale weak-signal areas, we select two decaying plage
regions that extend across the central meridian, and sample the B$\varphi$ on
either side (light/dark blue and green boxes in Figure 4). We assume that the
field vectors in these structures should not vary across the central meridian
per se. Figure 4 (right) shows the distribution of the zonal (East-West)
magnetic field in the area of negative (blue colors) and positive (green
colors) polarity flux. The distributions show a clear offset: for pixels
situated East of central meridian, the mean is positive, while for pixels in
the Western hemisphere, it is negative. Such an offset supports the notion
that a magnetic field situated in a region of decaying magnetic flux in the
Eastern hemisphere, on average, shows inclination patterns away from the solar
central meridian. Using observations with a lower noise level (5760 s) makes
the distributions narrower, but the offset is still present. This behavior is
in agreement with the explanation presented in Figure 2.

##  4 Effects of azimuth disambiguation and filling factor

For a simplified case of Bθ = 0, shown in Figure 2a, the transformation from
longitudinal and transverse components, measured in the image plane at a
heliographic longitude (central meridian distance) $\phi$, to heliographic
components Br and Bφ can be written as

$B_{r}=B_{||}\cos\phi\pm B_{\perp}\sin\phi=\frac{C_{||}}{\alpha}~{}{\rm
V}\cos\phi\pm\frac{C_{\perp}}{\sqrt{\alpha}}\sqrt[4]{{\rm Q}^{2}+{\rm
U}^{2}}\sin\phi$ (3)

$B_{\varphi}=-B_{||}\sin\phi\pm
B_{\perp}\cos\phi=-\frac{C_{||}}{\alpha}~{}{\rm
V}\sin\phi\pm\frac{C_{\perp}}{\sqrt{\alpha}}\sqrt[4]{{\rm Q}^{2}+{\rm
U}^{2}}\cos\phi.$ (4)

For this simple configuration, the $\pm$ ambiguity in Equations 3 and 4 is
resolved by requiring the two components to have the same sign in Equation 3
and consequently opposite signs in Equation 4. The first component (Equation
4) is positive East of the central meridian ($\phi<0$), and it is negative
West of the central meridian. The second component is always positive. For the
case when the field is mostly radial

$|\frac{-C_{||}~{}V}{\alpha}|>>\frac{C_{\perp}\sqrt[4]{Q^{2}+U^{2}}}{\sqrt{\alpha}}.$
(5)

Under this condition, the sum of the two components, which represents the sign
of Bφ, will be positive East of the central meridian and negative to the West.
However, in the case of a mostly horizontal field

$|\frac{-C_{||}~{}V}{\alpha}|<<\frac{C_{\perp}\sqrt[4]{Q^{2}+U^{2}}}{\sqrt{\alpha}}.$
(6)

The sign of Bφ will be always negative, independent of pixel location relative
to the central meridian. This is in qualitative agreement with the
observations mentioned in Section 1 that the disagreement in the sign of the
zonal (Bθ) and/or meridional (B$\varphi$) components of HMI and VSM
observations was only observed in pixels with weaker linear polarization
signals. In pixels with stronger linear polarization, the two instruments were
in agreement with each other with respect to the sign of Bφ and Bθ components.
Equations 5-6 suggest that the effect of sign reversal in Bφ could also depend
on magnetic filling factor $\alpha$. For some amplitudes of $|{-C_{||}~{}V}|$
and $C_{\perp}\sqrt[4]{Q^{2}+U^{2}}$ the inequity in Equation 5 could change
its sign if $\alpha$ is changed from being less then one to unity. The
appendix provides an example of a test done with SOLIS/VSM data, which
demonstrates that for $\alpha<1$, Bφ in three flux areas reverse their sign
when the area crosses the central meridian. When $\alpha$ is set to unity, the
same areas do not exhibit a sign reversal.

##  5 Discussion

We provide simple arguments and some observational evidence that a dichotomy
in properties of Zeeman polarization arising from longitudinal and transverse
field components, and our ability to interpret these signals especially in
unresolved structures (non-unity filling factor), combined with the azimuthal
disambiguation may lead to erroneous conclusions about the orientation of the
vector magnetic field, particularly in areas where the polarization signals
are weak. The systematic patterns may depend on the amplitude of noise and
whether the magnetic filling factor is resolved or assumed to be unity, and
thus, observations from different instruments could result in slightly
different orientation patterns (e.g., inclination of vector fields towards the
solar poles or towards the solar equator). Pixels with stronger polarization
signals are less affected (e.g., in sunspots, where $\alpha$ is relatively
large, and where Stokes Q, U, and V are typically strong), and thus, we expect
that only pixels with weaker polarization signals will show the effect
described above. Indeed, comparison of observations from SOLIS/VSM and SDO/HMI
show a good agreement between the two instruments in areas of strong fields
(e.g., sunspots), while in some areas of weak fields (e.g., plage), we see the
zonal (East-West) and/or meridional (North/South) components having opposite
sign. The amplitude of the effect will depend on the specific orientation of
the magnetic field vector relative to line-of-sight, and thus the impact of
this effect could vary across the solar disk in a somewhat complex way, making
it difficult (or perhaps impossible) to correct, since the signal inherently
originates from an unknown magnetic configuration. For an underlying radial-
field orientation, the effect could be mitigated by equalizing the amplitude
of the noise between the inferred longitudinal and transverse field
components. However, given the significant differences in noise levels between
inferred B|| and B⟂ such noise equalization could be done only at the expense
of increasing noise in the longitudinal field measurements (e.g., by using a
modulation schema allowing much longer integration time for states
corresponding to Stokes Q and U measurements). A test conducted using existing
HMI data shows that simply improving S/N for both transverse and longitudinal
fields does not completely eliminated the problem; while the amplitude of a
systematic inclination decreases slightly, the overall effect still remains.
We note that the systematic bias in the horizontal component of the magnetic
field may, in principle, lead to the effects similar in appearance to ones
that we discuss here. While our simplified example discussed in Section 4 uses
filling factor, it might be extremely difficult or even counterproductive to
model the exact magnetic and non-magnetic contribution in pixels outside
sunspots. The methods used to estimate the filling factor may inadvertently
introduce additional errors related to the unknown difference in line profiles
between magnetized and non-magnetized atmospheres in each pixel. The effects
of the contribution function in spatial direction (e.g., whether non-
magnetized component is predominantly located on one side of a pixel or it is
uniformly distributed around it as in Harker, 2017) are also unknown. Perhaps,
some of these issues may be addressed via spatially coupled inversion of
spectro-polarimetric image data as in van Noort (2012). That technique,
however, requires the knowledge of the Point-Spread-Function (PSF), which
could be achieved either by measuring the atmospheric seeing during the
observations (van Noort, 2017) or using adaptive optics. Thus, the method is
not applicable for the existing SOLIS/VSM observations.

Figure 4:  (Upper-left): half-disk image of radial Br field (white/black
correspond to positive/negative polarity), scaled to $\pm$200 G. Outlines show
pixels with the minimum-energy disambiguation applied (confid$\\_$disambig
keyword $\geq$ 60), which are included in the distributions shown in right
columns. (Lower-left: half-disk image of zonal Bφ, scaled to $\pm$100 G.
Region 1 (blue boxes) corresponds to area of negative polarity flux in
decaying flux region. Region 2 (green boxes) corresponds to a similar area but
of positive polarity flux. Boxes outlined by lighter/darker color are located
West/East of the central meridian. (Middle column): distribution of the zonal
(Bφ) component of the magnetic field in the Region 1 for two (720 s and 5760
s) averaging. (Right column): similar to middle column but for Region 2. Note:
salmon (light pinkish-orange) contours that outline pixels with the minimum-
energy disambiguation mask underlying black/white magnetic fields. To see both
the magnetic field and the contours, the PDF figure needs to be magnified by
300% or more.

###### Acknowledgements.

The authors thank Robert Cameron for his thoughtful comments, which helped
improving this article. We acknowledge help by Alexander Pevtsov in
preparation of the manuscript. This work is the result of discussions held at
three working meetings on the “Use of Vector Synoptic Maps for Modeling”.
Funding for these meetings was provided by the hosting groups at the
University of Oulu, Finland; the National Solar Observatory, USA; the Max
Planck Institute for Solar System Research, Germany; and by NASA’s Solar
Dynamics Observatory (Dr. Dean Pesnell). We also acknowledge the financial
support by the Academy of Finland to the ReSoLVE Centre of Excellence (project
no. 307411). AAP and LB acknowledge partial support by NASA grants
80NSSC17K0686 and NNX15AN43G. KDL acknowledges partial support from NASA/GSFC
grant 80NSSC19K0317. This work utilizes SOLIS data obtained by the NSO
Integrated Synoptic Program (NISP), managed by the National Solar Observatory.
HMI data used here are courtesy of NASA/SDO and the HMI science teams.

##  Appendix A The effect of magnetic fill factor on the horizontal components
of magnetic field

Figure 5 provides an example of a test done with SOLIS/VSM data, which
demonstrates that for $\alpha<1$, Bφ in three flux areas reverse its sign when
the area crosses central meridian. When $\alpha$ is set to unity, the same
areas do not exhibit the sign reversal.

Figure 5: Full disk images of Br (upper row), Bθ (middle) and Bφ (low row).
First and third columns show the vector field components in image coordinate
system, while second and fourth columns are in the heliographic coordinate
system. Three examples are marked by arrows of different color. Setting
$\alpha$ to unity emphasizes pixels, which otherwise are below the noise
threshold.

## References

  * Balasubramaniam and Pevtsov (2011) Balasubramaniam, K. S., and A. Pevtsov, 2011.  Ground-based synoptic instrumentation for solar observations.  In Proc. SPIE, vol. 8148 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , 814809.  10.1117/12.892824.
  * Borrero et al. (2011) Borrero, J. M., S. Tomczyk, M. Kubo, H. Socas-Navarro, J. Schou, S. Couvidat, and R. Bogart, 2011.  VFISV: Very Fast Inversion of the Stokes Vector for the Helioseismic and Magnetic Imager.  _Sol. Phys._ , 273(1), 267–293.  10.1007/s11207-010-9515-6, 0901.2702.
  * Cacciani et al. (1990) Cacciani, A., J. Varsik, and H. Zirin, 1990.  Observations of vector magnetic fields with a magneto-optic filter.  _Sol. Phys._ , 125, 173–178.  10.1007/BF00154786.
  * Centeno et al. (2014) Centeno, R., J. Schou, K. Hayashi, A. Norton, J. T. Hoeksema, Y. Liu, K. D. Leka, and G. Barnes, 2014.  The Helioseismic and Magnetic Imager (HMI) Vector Magnetic Field Pipeline: Optimization of the Spectral Line Inversion Code.  _Sol. Phys._ , 289(9), 3531–3547.  10.1007/s11207-014-0497-7, 1403.3677.
  * Couvidat et al. (2016) Couvidat, S., J. Schou, J. T. Hoeksema, R. S. Bogart, R. I. Bush, T. L. Duvall, Y. Liu, A. A. Norton, and P. H. Scherrer, 2016.  Observables Processing for the Helioseismic and Magnetic Imager Instrument on the Solar Dynamics Observatory.  _Sol. Phys._ , 291(7), 1887–1938.  10.1007/s11207-016-0957-3, 1606.02368.
  * Hagyard (1985) Hagyard, M. J., ed., 1985.  Measurements of Solar Vector Magnetic Fields, NASA Conference Proceedings 2375. Washington, DC: NASA.
  * Harker (2017) Harker, B. J., 2017.  VFISV Inversion Code Documentation for SOLIS/VSM Pipeline Implementation.  _Tech. rep._ , National Solar Observatory: Tucson, AZ.  1707.00002.
  * Hoeksema et al. (2014) Hoeksema, J. T., Y. Liu, K. Hayashi, X. Sun, J. Schou, et al., 2014.  The Helioseismic and Magnetic Imager (HMI) Vector Magnetic Field Pipeline: Overview and Performance.  _Sol. Phys._ , 289(9), 3483–3530.  10.1007/s11207-014-0516-8, 1404.1881.
  * Ichimoto et al. (1993) Ichimoto, K., T. Sakurai, E. Hiei, Y. Nishino, K. Shinoda, et al., 1993\.  Solar Flare Telescope project.  _Report of the National Astronomical Observatory of Japan_ , 1(4), 375–390.
  * Jefferies et al. (1989) Jefferies, J., B. W. Lites, and A. Skumanich, 1989.  Transfer of Line Radiation in a Magnetic Field.  _ApJ_ , 343, 920–935.
  * Jefferies and Mickey (1991) Jefferies, J., and D. L. Mickey, 1991.  On the Inference of Magnetic Field Vectors from Stokes Profiles.  _ApJ_ , 372, 694–702.
  * Keller et al. (2003) Keller, C. U., J. W. Harvey, and M. S. Giampapa, 2003.  SOLIS: an innovative suite of synoptic instruments.  In S. L. Keil and S. V. Avakyan, eds., Innovative Telescopes and Instrumentation for Solar Astrophysics, vol. 4853 of _Proc. SPIE_ , 194–204.  10.1117/12.460373.
  * Leka et al. (2009) Leka, K. D., G. Barnes, A. D. Crouch, T. R. Metcalf, G. A. Gary, J. Jing, and Y. Liu, 2009.  Resolving the 180∘ Ambiguity in Solar Vector Magnetic Field Data: Evaluating the Effects of Noise, Spatial Resolution, and Method Assumptions.  _Sol. Phys._ , 260(1), 83–108.  10.1007/s11207-009-9440-8.
  * Leka and Skumanich (1999) Leka, K. D., and A. Skumanich, 1999.  On the value of ‘$\alpha$AR’ from vector magnetograph data - I. Methods and Caveats.  _Sol. Phys._ , 188(1), 3–19.  10.1023/A:1005108632671.
  * Metcalf (1994) Metcalf, T. R., 1994.  Resolving the 180-degree ambiguity in vector magnetic field measurements: The ‘minimum’ energy solution.  _Sol. Phys._ , 155(2), 235–242.  10.1007/BF00680593.
  * Metcalf et al. (2006) Metcalf, T. R., K. D. Leka, G. Barnes, B. W. Lites, M. K. Georgoulis, et al., 2006.  An Overview of Existing Algorithms for Resolving the 180∘ Ambiguity in Vector Magnetic Fields: Quantitative Tests with Synthetic Data.  _Sol. Phys._ , 237(2), 267–296.  10.1007/s11207-006-0170-x.
  * Norton et al. (2006) Norton, A. A., J. P. Graham, R. K. Ulrich, J. Schou, S. Tomczyk, et al., 2006.  Spectral Line Selection for HMI: A Comparison of Fe I 6173 Å and Ni I 6768 Å.  _Sol. Phys._ , 239(1-2), 69–91.  10.1007/s11207-006-0279-y, astro-ph/0608124.
  * Pesnell et al. (2012) Pesnell, W. D., B. J. Thompson, and P. C. Chamberlin, 2012.  The Solar Dynamics Observatory (SDO).  _Sol. Phys._ , 275(1-2), 3–15.  10.1007/s11207-011-9841-3.
  * Pflug and Grigoryev (1986) Pflug, K., and V. M. Grigoryev, 1986.  The chain of the solar magnetographs and its results.  _Contributions of the Astronomical Observatory Skalnate Pleso_ , 15, 453–468.
  * Ronan et al. (1987) Ronan, R. S., D. L. Mickey, and F. Q. Orrall, 1987.  The derivation of vector magnetic fields from stokes profiles: Integral versus least squares fitting techniques.  _Sol. Phys._ , 113(1-2), 353–359.  10.1007/BF00147722.
  * Rudenko and Anfinogentov (2014) Rudenko, G. V., and S. A. Anfinogentov, 2014.  Very Fast and Accurate Azimuth Disambiguation of Vector Magnetograms.  _Sol. Phys._ , 289(5), 1499–1516.  10.1007/s11207-013-0437-y.
  * Scherrer et al. (2012) Scherrer, P. H., J. Schou, R. I. Bush, A. G. Kosovichev, R. S. Bogart, et al., 2012.  The Helioseismic and Magnetic Imager (HMI) Investigation for the Solar Dynamics Observatory (SDO).  _Sol. Phys._ , 275, 207–227.  10.1007/s11207-011-9834-2.
  * Severny (1965) Severny, A. B., 1965.  Introductory report.  In R. Lust, ed., Stellar and Solar Magnetic Fields, vol. 22 of _IAU Symposium_ , 238–266. North-Holland Pub. Co., Amsterdam.
  * Stenflo (2013) Stenflo, J. O., 2013.  Solar magnetic fields as revealed by Stokes polarimetry.  _A &A Rev._, 21, 66.  10.1007/s00159-013-0066-3, 1309.5454.
  * Stepanov and Severny (1962) Stepanov, V. E., and A. B. Severny, 1962.  A photoelectric method for measurements of the magnitude and direction of the solar magnetic field.  _Bull Crim Astrophys Obs_ , 28, 166–193.
  * van Noort (2012) van Noort, M., 2012.  Spatially coupled inversion of spectro-polarimetric image data. I. Method and first results.  _A &A_, 548, A5.  10.1051/0004-6361/201220220, 1210.4636.
  * van Noort (2017) van Noort, M., 2017.  Image restoration of solar spectra.  _A &A_, 608, A76.  10.1051/0004-6361/201731339, 1711.09629.
  * Virtanen et al. (2019) Virtanen, I. I., A. A. Pevtsov, and K. Mursula, 2019.  Structure and evolution of the photospheric magnetic field in 2010-2017: comparison of SOLIS/VSM vector field and BLOS potential field.  _Astronomy and Astrophysics_ , 624, A73.  10.1051/0004-6361/201834895, 1904.10740.


# Learning physics confers pose-sensitivity in structure-based virtual
screening

Pawel Gniewek
Atomwise Inc
San Francisco, CA, USA
<EMAIL_ADDRESS>
&Bradley Worley
Atomwise Inc
San Francisco, CA, USA
<EMAIL_ADDRESS>
&Kate Stafford
Atomwise Inc
San Francisco, CA, USA
<EMAIL_ADDRESS>
&Henry van den Bedem
Atomwise Inc
University of California, San Francisco
San Francisco, CA, USA
<EMAIL_ADDRESS>
&Brandon Anderson
Atomwise Inc
San Francisco, CA, USA
<EMAIL_ADDRESS>

###### Abstract

In drug discovery, structure-based virtual high-throughput screening (vHTS)
campaigns aim to identify bioactive ligands or “hits” for therapeutic protein
targets from docked poses at specific binding sites. However, while generally
successful at this task, many deep learning methods are known to be
insensitive to protein-ligand interactions, decreasing the reliability of hit
detection and hindering discovery at novel binding sites. Here, we overcome
this limitation by introducing a class of models with two key features: 1) we
condition bioactivity on pose quality score, and 2) we present poor poses of
true binders to the model as negative examples. The conditioning forces the
model to learn details of physical interactions. We evaluate these models on a
new benchmark designed to detect pose-sensitivity.

##  1 Introduction

Vast, make-on-demand chemical libraries like ENAMINE or Mcule have transformed
the scale of pharmaceutical, structure-based virtual high-throughput screening
(vHTS) campaigns [1]. To identify a “hit” from a library of candidate
molecules, structure-based vHTS methods predict binding affinity between a
protein and a ligand from their docked, bound complex, thereby assuming that
experimentally observed affinities correlate with protein-ligand interactions.
Conventional methods use empirical, physics-based approaches, which attempt to
calculate the binding free energy of complex formation. By contrast newer deep
learning (DL) approaches are trained on large data sets using implicit
features to predict activity. These statistical models can outperform physics-
based approaches in retrospective tests for predicting activity.

Early structure-based DL methods in vHTS centered on Convolutional Neural
Networks (CNNs), representing protein-ligand structures by a 3D voxel grid to
predict activity [2, 3, 4, 5]. Although generally effective [6], a drawback of
CNNs is that they are not rotationally invariant, and require more parameters
than alternative representations. Consequently, Graph Convolutional Networks
(GCNs) [7, 8], and their position-based generalizations [9, 10, 11, 12, 13,
14, 15, 16], have gained popularity by capitalizing on physical symmetries and
using data more efficiently.

Recent studies have suggested that the performance of structure-based machine
learning methods is partly driven by proteochemometric (PCM)-like features [5,
17, 18]. Rather than responding to specific interactions between the ligand
and the binding site, the model learns independent ligand and protein
signatures. This deficiency manifests as a drop in predictive performance when
the model is confronted with a previously unseen binding site on the same
protein, especially when that site partially overlaps with a canonical site
and portions of the canonical site lie within the receptive field of the
network. Thus the model can be biased by detected but chemically irrelevant
features provided in the input data – analogous to the “Picasso Problem” in 2D
CNNs. For example, the model may highly rank ATP-competitive binders for a
neighboring allosteric site on a protein kinase. This limitation critically
hinders discovery of new chemical matter, or the ability to target novel sites
on proteins.

Simultaneous training on ligand pose quality and affinity can improve pose
sensitivity [19]. Here, we build on that observation and present a multi-task
architecture which simultaneously evaluates bioactivity and the physics-based
Vina [20] score of the pose. We then condition its bioactivity task on a
separate measure of pose quality, as defined by the score assigned to a pose
by PoseRanker, a GCN developed to identify high-quality poses from docking
[13]. Finally, we use poor poses of true positive compounds as additional
negative examples. We demonstrate that our architecture improves pose-
sensitivity on several benchmarks.

##  2 Methods

###  2.1 Neural Network Architectures

Our model is a GCN-based architecture with atoms as vertices, and pair-wise,
distance dependent edges. We consider only receptor atoms within 7Å of any
ligand atom. In the first two graph convolutional layers any two atoms share
an edge if they are within 4Å of each other. We then extract ligand-only
features and follow with two ligand-only layers. We sum-pool the final ligand-
only layer to create an embedding for the multi-task multilayer perceptrons at
the top of the network. The network architecture is further described in
Sections A.1 and A.2. We predict three outputs using these embeddings:
bioactivity, the PoseRanker pose quality score, and the Vina docking score.
This is performed in two stages. First, we compute the PoseRanker and Vina
score prediction by passing the embedding through two independent multilayer
perceptrons. We then form a conditioned embedding by concatenating the input
embedding with the PoseRanker score prediction, which is passed to a third
multilayer perceptron to compute the activity prediction [21]. Section A.3
describes model training parameters.

Figure 1: Multi-task architecture of the MT-4b model. Table 1: Deep Learning Neural Networks used in this work. Architectures are either a 3D CNN-based model or a GCN-based model. Multi-task (MT) models include additional Vina [20, 22] or PoseRanker [13] tasks, and may include additional conditioning layer.  Name | Architecture | Vina | PoseRanker | Pose-negatives | Cond.
---|---|---|---|---|---
CNN | 3D CNN based | ✗ | ✗ | ✗ | ✗
GCN | GCN based | ✗ | ✗ | ✗ | ✗
MT-1 | GCN based | ✓ | ✗ | ✗ | ✗
MT-2 | GCN based | ✓ | ✓ | ✗ | ✗
MT-3 | GCN based | ✓ | ✓ | ✗ | ✓
MT-4a | GCN based | ✓ | ✓ | ✓ | ✗
MT-4b | GCN based | ✓ | ✓ | ✓ | ✓

###  2.2 Data

The training dataset consists of binding affinity measurements collected from
publicly available sources (Chembl, Pubchem) and commercial databases (Reaxys,
Liceptor). We considered quantitative measurements only with
$\textrm{pKi}\in(0,11)$. We labeled compounds as active if their measured pKi
(or pIC50) is less than 10$\mu$M; otherwise they were labeled as inactive.
There are more active than inactive compounds, so to balance the training set
we randomly assigned each of the active compounds as inactive decoys for
another, dissimilar protein target. For some models we also used pose-
negatives - poor poses of active compounds treated as negatives (Section A.4).
We selected a set of 12 diverse proteins (D12) to serve as a hold-out test set
and excluded all close homologs of the D12 proteins (>95% sequence similarity)
from training. Training set covers more than 3800 diverse proteins, and counts
4.8M (5.8M) data points without (with) pose-negatives. The hold-out set counts
about 33,000 compounds and 12 proteins. We docked compounds with the GPU-
enabled docking engine, CUina [22]. We used the PoseRanker [13] top-ranked
pose for scoring.

###  2.3 Numerical experiments

To evaluate pose-sensitivity of our models, we scored each of the active
target-compound pairs in D12 three times using: i) the top pose, ii) a poor
pose, iii) the average score of four physically implausible poses, each
obtained by randomly rotating the top pose around its center of mass. (See
Figure 2a.) Our measure of pose-sensitivity is the median of the drop of the
bioactivity score between good and poor/implausible poses. We expect that a
model sensitive to ligand-receptor interactions reveals sharply reduced
bioactivity scores for poor or implausible poses compared to good poses.

Next, we designed a benchmark based on a protein with multiple binding sites
to investigate the Picasso Problem described above. To monitor how the
neighbor binding sites interfere with the model’s inference, we selected three
binding sites on the same protein kinase hZAP70: the ATP site, a nearby
allosteric site 6-10Å away, and a distant SH2 site over 50Å away (Figure 2).
We selected a diverse set of known ATP-site-binding kinase inhibitors (about
300 compounds labeled as actives) and mixed them with $10^{5}$ randomly
selected compounds from a screening library (MCULE (17/08/2018); labeled as
inactive). Each compound was docked to each site. We expect that a good model
can distinguish kinase inhibitors against a background of random molecules
docked to the ATP binding site (ROC AUC $\gg$ 0.5). However, a pose-sensitive
model should not be biased by the neighbor ATP-site when compounds are docked
to the allosteric site, and therefore should have near random performance. To
account for possible biases in the training set, we also examined the distant
control site on the SH2 domain, which likewise does not bind ATP-site targeted
molecules and is expected to have near-random performance in identifying ATP-
site binders.

Figure 2: (a) A schematic of the pose-sensitivity experiment. For each target-
ligand pair we generate three types of poses: i) a good pose ii) a poor pose
and iii) an implausible pose. Good and poor poses are the highest and lowest
ranked poses by PoseRanker after sampling with CUina. Implausible poses are
obtained from good poses by a random rotation of the ligand around its center
of mass. Scores for the poor and implausible poses are subtracted from the
score of the good pose. (b) Human ZAP 70 protein (PDB: 2ozo) indicating ATP
site (grey), allosteric site (red), and spatially distant SH2 site (blue). A
set of $10^{5}$ compounds is docked to each site and scored with the Deep
Learning models.

##  3 Results

Figure 3: (a) Median probability drops between good and poor poses. (b) Median
probability drops between good and implausible poses. (c) Picasso problem
experiment. $10^{5}$ diverse compounds (labeled as 0, non-binders) mixed with
c.a. 300 kinase inhibitors (labeled as 1, binders) were docked and scored with
at three binding sites i) ATP binding site ii) allosteric binding site iii)
binding site at the SH2 domain, See Figure 2b.

Figure S3 shows that our baseline models perform well on the holdout set, with
GCN-based models slightly outperforming CNN. However, both single-task models
enrich ATP-site compounds at the allosteric site (Figure 3). Because the ATP
site is in the receptive field of these two networks when evaluating the
allosteric site, GCN and CNN models make predictions as if they were shown the
ATP site (Figure 2b). This cannot be due to a biased training set, as
screening the spatially distant SH2 site gave no enrichment of ATP-site
compounds (Figure 3c). In other words, the single-task CNN and GCN models are
effectively PCMs. Further, they are insensitive not only to poor poses, but
also to implausible poses (Figure 3a,b). This behavior was observed previously
with 3D-grid-based CNNs [4, 19], but a general solution has been lacking.

At first it is unintuitive that the single task activity models trained on
structural data do not use that structural information about ligand-receptor
interactions. Apparently, the model has no direct incentive to learn the
structural features of the ligand-receptor interactions, so models often
neglect them. Training a multi-task model where the additional tasks require
embeddings sensitive to structure should alleviate the problem. Indeed, the
Multi-Task (MT) models show increased pose-sensitivity. Adding a docking score
regression task (model MT-1) leads to a model that penalizes implausible poses
and decreases the enrichment of kinase inhibitors at the allosteric site of
the hZAP70 protein (Figure 3b,c). However, because poor poses are generally
far from native but do not have any atomic clashes, the MT-1 model is unable
to distinguish between good and poor poses (Figure 3a). Interestingly, adding
pose quality regression as the third task (MT-2) or conditioning the activity
task on the pose quality (MT-3) do not increase sensitivity to poor poses.
Models shown only top-ranked poses do not learn what bad poses look like. To
correct for this, we introduced "pose negatives": bottom-ranked poses of true
actives treated as negative examples. The resulting models (MT-4a and MT-4b)
penalize physically implausible poses as well as poor poses. The same models
also mitigate the Picasso problem. The model that does not condition its
activity task on pose quality (MT-4a) is still somewhat more prone to the
Picasso problem compared to the conditioning model (MT-4b).

##  4 Conclusions

Multi-task architectures can lead to models that are capable of predicting
bioactivity of compounds while making full use of the structural data provided
for the inference. Forcing the model to learn orthogonal tasks regularizes the
final model. Pose sensitive models enable structure-based deep learning
approaches in drug discovery to target novel binding sites, and unlock
previously undruggable proteins.

## Acknowledgments and Disclosure of Funding

We thank Michael Mysinger and Adrian Stecula for helpful comments.

## References

  * Irwin and Shoichet [2016] John J. Irwin and Brian K. Shoichet.  Docking Screens for Novel Ligands Conferring New Biology: Miniperspective.  _Journal of Medicinal Chemistry_ , 59(9):4103–4120, May 2016.  ISSN 0022-2623, 1520-4804.  doi: 10.1021/acs.jmedchem.5b02008.  URL https://pubs.acs.org/doi/10.1021/acs.jmedchem.5b02008.
  * Wallach et al. [2015] Izhar Wallach, Michael Dzamba, and Abraham Heifets.  AtomNet: A Deep Convolutional Neural Network for Bioactivity Prediction in Structure-based Drug Discovery.  _arXiv:1510.02855 [cs, q-bio, stat]_ , October 2015.  URL http://arxiv.org/abs/1510.02855.  arXiv: 1510.02855.
  * Ragoza et al. [2017] Matthew Ragoza, Joshua Hochuli, Elisa Idrobo, Jocelyn Sunseri, and David Ryan Koes.  Protein–Ligand Scoring with Convolutional Neural Networks.  _Journal of Chemical Information and Modeling_ , 57(4):942–957, April 2017.  ISSN 1549-9596, 1549-960X.  doi: 10.1021/acs.jcim.6b00740.  URL https://pubs.acs.org/doi/10.1021/acs.jcim.6b00740.
  * Stepniewska-Dziubinska et al. [2018] Marta M Stepniewska-Dziubinska, Piotr Zielenkiewicz, and Pawel Siedlecki.  Development and evaluation of a deep learning model for protein–ligand binding affinity prediction.  _Bioinformatics_ , 34(21):3666–3674, November 2018.  ISSN 1367-4803, 1460-2059.  doi: 10.1093/bioinformatics/bty374.  URL https://academic.oup.com/bioinformatics/article/34/21/3666/4994792.
  * Boyles et al. [2019] Fergus Boyles, Charlotte M Deane, and Garrett M Morris.  Learning from the ligand: using ligand-based features to improve binding affinity prediction.  _Bioinformatics_ , page btz665, August 2019.  ISSN 1367-4803, 1460-2059.  doi: 10.1093/bioinformatics/btz665.  URL https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btz665/5554651.
  * Hsieh et al. [2019] Chung-Han Hsieh, Li Li, Roeland Vanhauwaert, Kong T. Nguyen, Mary D. Davis, Guojun Bu, Zbigniew K. Wszolek, and Xinnan Wang.  Miro1 Marks Parkinson’s Disease Subset and Miro1 Reducer Rescues Neuron Loss in Parkinson’s Models.  _Cell Metabolism_ , 30(6):1131–1140.e7, December 2019.  ISSN 15504131.  doi: 10.1016/j.cmet.2019.08.023.  URL https://linkinghub.elsevier.com/retrieve/pii/S1550413119304978.
  * Kipf and Welling [2017] Thomas N. Kipf and Max Welling.  Semi-Supervised Classification with Graph Convolutional Networks.  _arXiv:1609.02907 [cs, stat]_ , February 2017.  URL http://arxiv.org/abs/1609.02907.  arXiv: 1609.02907.
  * Gilmer et al. [2017] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl.  Neural message passing for Quantum chemistry.  In _Proceedings of the 34th International Conference on Machine Learning - Volume 70_ , ICML’17, pages 1263–1272, Sydney, NSW, Australia, August 2017. JMLR.org.
  * Behler and Parrinello [2007] Jörg Behler and Michele Parrinello.  Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces.  _Physical Review Letters_ , 98(14):146401, April 2007.  doi: 10.1103/PhysRevLett.98.146401.  URL https://link.aps.org/doi/10.1103/PhysRevLett.98.146401.
  * Schütt et al. [2017] Kristof T. Schütt, Pieter-Jan Kindermans, Huziel E. Sauceda, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert Müller.  SchNet: A continuous-filter convolutional neural network for modeling quantum interactions.  _arXiv:1706.08566 [physics, stat]_ , December 2017.  URL http://arxiv.org/abs/1706.08566.  arXiv: 1706.08566.
  * Feinberg et al. [2018] Evan N. Feinberg, Debnil Sur, Zhenqin Wu, Brooke E. Husic, Huanghao Mai, Yang Li, Saisai Sun, Jianyi Yang, Bharath Ramsundar, and Vijay S. Pande.  PotentialNet for Molecular Property Prediction.  _ACS Central Science_ , 4(11):1520–1530, November 2018.  ISSN 2374-7943, 2374-7951.  doi: 10.1021/acscentsci.8b00507.  URL https://pubs.acs.org/doi/10.1021/acscentsci.8b00507.
  * Lim et al. [2019] Jaechang Lim, Seongok Ryu, Kyubyong Park, Yo Joong Choe, Jiyeon Ham, and Woo Youn Kim.  Predicting Drug–Target Interaction Using a Novel Graph Neural Network with 3D Structure-Embedded Graph Representation.  _Journal of Chemical Information and Modeling_ , 59(9):3981–3988, September 2019.  ISSN 1549-9596, 1549-960X.  doi: 10.1021/acs.jcim.9b00387.  URL https://pubs.acs.org/doi/10.1021/acs.jcim.9b00387.
  * Stafford et al. [2021] Kate Stafford, Brandon M. Anderson, Jon Sorenson, and Henry van den Bedem.  AtomNet PoseRanker: Enriching Ligand Pose Quality for Dynamic Proteins in Virtual High Throughput Screens.  September 2021.  doi: 10.33774/chemrxiv-2021-t6xkj.  URL https://chemrxiv.org/engage/chemrxiv/article-details/614b905e39ef6a1c36268003.
  * Thomas et al. [2018] Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley.  Tensor field networks: Rotation- and translation-equivariant neural networks for 3D point clouds.  _arXiv:1802.08219 [cs]_ , May 2018.  URL http://arxiv.org/abs/1802.08219.  arXiv: 1802.08219.
  * Anderson et al. [2019] Brandon Anderson, Truong Son Hy, and Risi Kondor.  Cormorant: Covariant Molecular Neural Networks.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://papers.nips.cc/paper/2019/hash/03573b32b2746e6e8ca98b9123f2249b-Abstract.html.
  * Townshend et al. [2021] Raphael J. L. Townshend, Martin Vögele, Patricia Suriana, Alexander Derry, Alexander Powers, Yianni Laloudakis, Sidhika Balachandar, Bowen Jing, Brandon Anderson, Stephan Eismann, Risi Kondor, Russ B. Altman, and Ron O. Dror.  ATOM3D: Tasks On Molecules in Three Dimensions.  _arXiv:2012.04035 [physics, q-bio]_ , June 2021.  URL http://arxiv.org/abs/2012.04035.  arXiv: 2012.04035.
  * Sieg et al. [2019] Jochen Sieg, Florian Flachsenberg, and Matthias Rarey.  In Need of Bias Control: Evaluating Chemical Data for Machine Learning in Structure-Based Virtual Screening.  _Journal of Chemical Information and Modeling_ , 59(3):947–961, March 2019.  ISSN 1549-9596, 1549-960X.  doi: 10.1021/acs.jcim.8b00712.  URL https://pubs.acs.org/doi/10.1021/acs.jcim.8b00712.
  * Chen et al. [2019] Lieyang Chen, Anthony Cruz, Steven Ramsey, Callum J. Dickson, Jose S. Duca, Viktor Hornak, David R. Koes, and Tom Kurtzman.  Hidden bias in the DUD-E dataset leads to misleading performance of deep learning in structure-based virtual screening.  _PLOS ONE_ , 14(8):e0220113, August 2019.  ISSN 1932-6203.  doi: 10.1371/journal.pone.0220113.  URL https://dx.plos.org/10.1371/journal.pone.0220113.
  * Francoeur et al. [2020] Paul G. Francoeur, Tomohide Masuda, Jocelyn Sunseri, Andrew Jia, Richard B. Iovanisci, Ian Snyder, and David R. Koes.  Three-Dimensional Convolutional Neural Networks and a Cross-Docked Data Set for Structure-Based Drug Design.  _Journal of Chemical Information and Modeling_ , 60(9):4200–4215, September 2020.  ISSN 1549-9596, 1549-960X.  doi: 10.1021/acs.jcim.0c00411.  URL https://pubs.acs.org/doi/10.1021/acs.jcim.0c00411.
  * Trott and Olson [2010] Oleg Trott and Arthur J. Olson.  AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading.  _Journal of Computational Chemistry_ , 31(2):455–461, 2010.  ISSN 1096-987X.  doi: https://doi.org/10.1002/jcc.21334.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.21334.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.21334.
  * Long et al. [2018] Mingsheng Long, Zhangjie Cao, Jianmin Wang, and Michael I. Jordan.  Conditional Adversarial Domain Adaptation.  _arXiv:1705.10667 [cs]_ , December 2018.  URL http://arxiv.org/abs/1705.10667.  arXiv: 1705.10667.
  * Morrison et al. [2020] Adrian Morrison, Greg Friedland, and Izhar Wallach.  CUina: An Efficient GPU Implementation of AutoDock Vina.  August 2020.  URL https://blog.atomwise.com/efficient-gpu-implementation-of-autodock-vina.

##  Appendix A Appendix

###  A.1 Graph Neural Network architecture

Throughout the paper we use a graph convolutional network (GCN) to generate
the embeddings used in our conditional multi-tasking architecture. Our GCN
takes the three dimensional coordinates of the protein-ligand pose, along with
a one-hot atom encoding that simultaneously identifies element type,
protein/ligand, and hybridization state. We define our connectivity purely by
radial functions without use of chemical bonds. See Figure S1 for details.

We use a radial cutoff $R^{l}_{c}$ at layer $l$ to define a radial graph, with
the neighborhood for atom $i$ defined as
$\mathcal{N}(i)=\\{j:d_{ij}<R_{c}\\}$, where
$d_{ij}=\|\mathbf{r}_{i}-\mathbf{r}_{j}\|$ is the pairwise distance between
atoms $i$ and $j$. At each layer $l$ we have a feature vector
$f^{l}_{i}\in\mathcal{R}^{N_{f}^{l}}$ at atom $i\in\mathcal{S}_{P/L/P+L}$,
where $N_{f}^{l}$ is the number of features at layer $l$, and
$\mathcal{S}_{P}$, $\mathcal{S}_{L}$,
$\mathcal{S}_{P+L}=\mathcal{S}_{P}\cup\mathcal{S}_{L}$ are the set of atoms in
the protein, ligand, or protein-ligand complex respectively. Our network can
be configured so that $P$, $L$, or $P+L$ atoms can be used as either source
atoms ($i$) or target atoms ($j$) on a layer-by-layer basis.

Our graph convolutional block is based upon a continuous filter convolution
[10], $f^{\rm conv}_{i}\leftarrow\sum_{j}W^{l}(d_{ij})f^{l}_{j}$. We use
convolutional kernels
$W^{l}(r)=\sum_{n}^{N_{b}^{l}}w_{n}^{l}j_{0}(rz_{0n}/R^{l}_{c})$ constructed
from a linear combinations of zeroth-order spherical Bessel functions
$j_{0}(x)=\sin x/x$. Here $z_{0n}$ is the $n$-th zero of $j_{0}(x)$, and
$w_{n}^{l}$ are learnable weights, and the number of basis functions is chosen
as $N_{b}^{l}=\lceil R_{c}^{l}/\delta_{l}\rceil$, where
$\delta_{l}=0.15625\textup{\AA}$. This parametrization ensures that the
convolutional kernels vanish at the edge of an internal neighborhood. Since we
are not calculating forces, we do not need a smooth cutoff here. After the
graph convolution step, we include a ${\rm Linear}-{\rm LeakyReLU}-{\rm
Linear}$ layer. We do not apply a residual connection at the output of a
convolution block. Instead, each convolutional layer takes as input all layers
using a bottleneck layer. We found this performed better empirically than
standard skip connections from the previous layer of the same set of source
atoms.

Our network has five total graph convolutional blocks. In the first two graph
conv blocks, we consider all ligand and receptor atoms so $i,j\in S_{P+L}$,
and choose $R_{c}=5\textup{\AA}$ and $N^{l}_{F}=64$ filters. The third graph
conv block, the cutoff radius and filters are increased to
$R_{c}=7\textup{\AA}$ and $N^{l}_{F}=128$, however only atoms in the ligand
are considered as destination atoms (i. e., $i\in S_{P+L}$, but $j\in S_{L}$
only). Finally, in the last two layers only ligand atom embeddings are
aggregated ($i,j\in S_{L}$) with a cutoff $R_{c}=7\textup{\AA}$ and
$N^{l}_{F}=128$.

To construct our final embedding, we first apply a readout operation to the
protein and ligand features independently at each layer: $z_{P}^{l}=\sum_{i\in
S_{P}}f_{i}^{l}$, $z_{L}^{l}=\sum_{i\in S_{L}}f_{i}^{l}$. We then concatenate
these embeddings together, ${\bf z}_{\rm
read}=\bigoplus_{l=0}^{1}z_{P}^{l}\oplus\bigoplus_{l=0}^{5}z_{L}^{l}$, and
finally construct ${\bf z}$ from two iterations of ${\rm Dropout}_{0.2}-{\rm
LeakyReLU}-{\rm Linear}$ applied to ${\bf z}_{\rm read}$. This final embedding
is then used as the input to the conditional multi-task prediction layer.

Figure S1: Architecture of the Graph Cronvolutional Network.

###  A.2 Conditional Multi-task architectures

The techniques used to condition network outputs on others (Figure 1) largely
follow those outlined by Long et al. [21]. Given the shared embedding
$\mathbf{z}_{\rm read}=g({f^{0...L}})$, the PoseRanker score
$\hat{y}_{\mathrm{pose}}$ and the Vina score $\hat{y}_{\mathrm{vina}}$ were
computed by passing $\mathbf{z}_{\rm read}$ through separate MLPs
$y_{\mathrm{pose}}$ and $y_{\mathrm{vina}}$, respectively. A conditioned
embedding $\mathbf{z}^{\prime}$ was then formed by passing $\mathbf{z}$
through the conditional map:

$h:\mathbf{z}_{\rm read}\mapsto\begin{pmatrix}\mathbf{z}_{\rm
read}\,\sigma(\hat{y}_{\mathrm{pose}})\\\ \mathbf{z}_{\rm
read}\left(1-\sigma(\hat{y}_{\mathrm{pose}})\right)\end{pmatrix}$

where $\sigma:x\mapsto(1+e^{-x})^{-1}$, and passed to a final MLP
$p_{\mathrm{active}}$ to compute the activity score. Alternatively,
$\mathbf{z}^{\prime}$ may be passed through a second conditional map,

$h^{\prime}:\mathbf{z}^{\prime}\mapsto\begin{pmatrix}\mathbf{z}^{\prime}\\\
\hat{y}_{\mathrm{vina}}\end{pmatrix}$

to obtain an embedding $\mathbf{z}^{\prime\prime}$ that has been conditioned
on both the PoseRanker score and the Vina score (Figure S2).

Figure S2: Active task conditioned on PoseRanker and Vina scores.

###  A.3 Training

For every architecture, we trained 6 models for 10 epochs, each with
$\nicefrac{{5}}{{6}}$th of the data as the training set, and
$\nicefrac{{1}}{{6}}$th left out for cross-fold validation. Each data cross-
fold contains clusters of proteins that share more than 70% of the sequence
similarity. We used the ADAM optimizer with the learning rate
$\mathrm{lr}=0.001$, and sampled targets with replacement, proportionally to
the number of active compounds associated with that target. Targets without
any measured active compounds were pruned from the training set.

###  A.4 Data augmentation with pose-negatives

We generated 64 ligand poses for every ligand-target pair with the CUina
docking engine. Next, we reranked poses with PoseRanker according to their
similarity to a native pose [13], and selected the top 16 poses. We used the
top-ranked pose as a ’good’ pose in training and scoring, whereas we used the
bottom-ranked (16th) pose as a ’poor’ (pose-negative) pose and considered
inactive (non-binder).

Figure S3:  ROC AUC performance on the D12 benchmark.