Cdr Clock Data Recovery Cadence

the data-parallel setting, where data is partitioned across different samples, whereas, CD is used in the model-parallelism setting, where data is partitioned across the parameter space. At the core of our solutions to both these algorithms is a method for Byzantine-resilient matrix-vector (MV) multiplication; and for that, we propose a method based on data encoding and error correction over real numbers to combat adversarial attacks.

⌋Corrupt Worker

nodes, which is information-theoretically optimal. We give deterministic guarantees, and our method does not assume any probability distribution on the data. We develop a sparse encoding scheme which enables computationally efficient data encoding and decoding.

We Demonstrate A Trade-Offbetween

the corruption threshold and the resource requirements (storage, computational, and communication

Complexity). As An Example, For T ≤M

3 , our scheme incurs only a constant overhead on these resources, over that required by the plain distributed PGD/CD algorithms which provide no adversarial protection. To the best of our knowledge, ours is the first paper that connects MV multiplication with CD and designs a specific encoding matrix for MV multiplication whose structure we can leverage to make CD secure against adversarial attacks.

Our encoding scheme extends efficiently to (i) the data streaming model, in which data samples come in an online fashion and are encoded as they arrive, and (ii) making stochastic gradient descent (SGD) Byzantine-resilient. In the end, we give experimental results to show the efficacy of our proposed schemes.

Ntroduction

Map-reduce architecture [DG08] is implemented in many distributed learning tasks, where there is one designated machine (called the master) that computes the model iteratively, based on the inputs from the worker machines at each iteration, typically using descent techniques, like (proximal) gradient descent, coordinate descent, stochastic gradient descent, the Newton’s method, etc. The worker nodes perform the required computations using local data, distributed to the nodes [ZWLS10]. Several other architectures, including having no hierarchy among the nodes have been explored [LZZ+17].

Arxiv:1907.02664V2 [Cs.Dc] 4 Nov 2020

In several applications of distributed learning, including the Internet of Battlefield Things (IoBT) [A+18], federated optimization [Kon17], the recruited worker nodes might be partially trusted with their computa- tion. Therefore, an important question is whether we can reliably perform distributed computation, taking advantage of partially trusted worker nodes. These Byzantine adversaries can collaborate and arbitrarily saries has a long history [LSP82], and there has been recent interest in applying this computational model to large-scale distributed learning [BMGS17,CWCP18,CSX17].

In this paper, we study Byzantine-tolerant distributed optimization to learn a regularized generalized linear model (GLM) (e.g., linear/ridge regression, logistic regression, Lasso, SVM dual, constrained mini- mization, etc.). We consider two frameworks for distributed optimization: (i) data-parallelism architecture, where data points are distributed across different worker nodes, and in each iteration, they all parallelly com- pute gradients on their local data and master aggregates them to update the parameter vector using gradient descent (GD) [BT89,Bot10,DCM+12]; and (ii) model-parallelism architecture, where data points are parti- tioned across features, and several worker nodes work in parallel, updating different subsets of coordinates of the model/parameter vector through coordinate descent (CD) [BKBG11,Wri15,RT16]. Note that GD requires full gradients to update the parameter vector; and if full gradients are too costly to compute, we can reduce the per-iteration cost by using CD,1 which also has been shown to be very effective for solving generalized lin- ear models, and is particularly widely used for sparse logistic regression, SVM, and Lasso [BKBG11]. Given its simplicity and effectiveness, CD can be chosen over GD in such applications [Nes12]. Computing gradients in the presence of Byzantine adversaries has been recently studied [BMGS17, CSX17, CWCP18, YCRB18, AAL18,SX19,XKG19,YCRB19,GV19,RWCP19,LXC+19,GHYR19,YLR+19,DD20b,DD20a,HKJ20], and we discuss them in detail Section 3 where we also put our work in context. However, as far as we know, making CD robust to Byzantine adversaries has not received much attention, and to the best of our knowl- edge, ours is the first paper that studies CD against Byzantine attacks and provides an efficient solution for that.

Our Contributions

We propose Byzantine-resilient distributed optimization algorithms both for PGD and CD based on data encoding and error correction (over real numbers). As mentioned above, there have been several papers that provide different methods for gradient computation in the presence of Byzantine adversaries, however, our proposed algorithm differs from them in one or more of the following aspects: (i) it does not make statistical assumptions on the data or Byzantine attack patterns; (ii) it can tolerate up to a constant fraction (< 1/2) of the worker nodes being Byzantine, which is information-theoretically optimal; and (iii) it enables a trade-off (in terms of storage and computation/communication overhead at the master and the worker nodes) with Byzantine adversary tolerance, without compromising the efficiency at the master node. We give the same guarantees for CD also.

First we design a coding scheme for distributed matrix-vector (MV) multiplication, specifically, for op- erating in the presence of Byzantine adversaries, and use that in both our algorithms for PGD and CD to and has been known for some time (see, for example, [LLP+18,DCG16]), however, it is not clear whether we can use MV multiplication methods for CD also. Indeed, since each CD update has a different requirement than that of gradient computation, a general-purpose algorithm for MV multiplication may not be applicable for CD. One distinction is that in gradient computation, we only need to encode the data to compute the MV multiplication, whereas, in CD, in addition to data encoding, since workers update few coordinates of different parts of the parameter vector in parallel, we need to encode the parameter vector as well for master to be able to decode that. In this paper, we design our encoding matrix for MV multiplication in such a way that it is sparse and has a regular structure of non-zero entries (see (11) for the encoding matrix for any worker), which makes it applicable for CD too. This leads to efficient solutions for both PGD and CD, 1Alternatively, we can also use SGD to reduce the per-iteration cost, and we give a method for making SGD Byzantine- resilient in Section 6.1.

which are our main focus in this paper. Inspired from the real-error correction (or sparse reconstruction) problem [CT05], we develop efficient encoding/decoding procedures for MV multiplication, where we encode the data matrix and distribute it to the m worker nodes, and to recover the MV product at the master, we reduce the decoding problem to the sparse reconstruction or real-error correction problem [CT05]. Note that in PGD, we only need to encode the data, whereas, in CD, we also need to encode the parameter vector, and our coding scheme should facilitate the requirement that the update on a small fraction of the encoded parameter vector should affect only a small fraction of the original parameter vector. This is a non-trivial requirement, and our coding scheme for MV multiplication is designed in such a way that it supports this requirement in an efficient manner; see Section 2.2 for a description on plain distributed CD, Section 2.5 for our approach to making CD robust to Byzantine attacks, and Section 5 for a complete solution for Byzantine-resilient CD. In the context of PGD/CD, for decoding, the master node processes the inputs from the worker nodes, either to compute the true gradient in the case of PGD or to facilitate the computation at the worker nodes in the summarized in Theorem 1 (on page 9) for PGD and Theorem 2 (on page 11) for CD, and demonstrate a trade-offbetween the Byzantine resilience (in terms of the number of adversarial nodes) and the resource requirement (storage, computational, and communication complexity).

, Our

scheme incurs only a constant overhead on these resources, over that required by the plain distributed PGD and CD algorithms which provide no adversarial protection. Our coding schemes can handle both Byzantine attacks and missing updates (e.g., caused by delay or asynchrony of worker nodes). Our encoding process is also efficient. Though data encoding is a one-time process, it has to be efficient to harness the advantage of distributed computation. We design a sparse encoding process, based on real-error correction, which enables efficient encoding, and the worker nodes encode data using the sparse structure. This allows encoding with

M

m−2t (which is a constant, even if t is a constant (< 1

) Fraction Of M), And A One-Time

total computation cost for encoding is O((1 + 2t)nd). Note that the time for data encoding is a factor of (1 + 2t) (where t is the corruption threshold) more than the time required for plain data distribution which is O(nd), the size of the data matrix.

We extend our encoding scheme in a couple of important ways: first, to make the stochastic gradient descent (SGD) algorithm Byzantine-resilient without compromising much on the resource requirements; and second, to handle streaming data efficiently, where data points arrives one by one (and we encode them as they arrive), rather than being available at the beginning of the computation; we also give few more applications of our method. For the streaming model, more specifically, our encoding requires the same amount of time, irrespective of whether we encode all the data at once, or we get data points one by one (or in batches) and we encode them as they arrive. This setting encompasses a more realistic scenario, in which we design our coding scheme with the initial set of data points and distribute the encoded data among the workers. Later on, when we get some more samples, we can easily incorporate them into our existing encoded setup. See Section 6 for details on these extensions.

Paper Organization

We present our problem formulation, description of the plain distributed PGD and CD algorithms, and the high-level ideas of our Byzantine-resilient algorithms for both PGD and CD along-with our main results in Section 2.

We give detailed related work in Section 3.

We Present Our Full Coding Schemes For Mv

multiplication and also for gradient computation for PGD along-with a complete analysis of their resource requirements in Section 4. In Section 5, we provide a complete solution to CD. In Section 6, we show how our method can be extended to SGD and to the data streaming model. We also discuss applicability of our method to a few more important applications in that section. In Section 7, we show numerical results of our method: we show the efficiency of our method for both gradient descent (GD) and coordinate descent (CD) by running them to solve linear regression on two datasets (moderate and large) and plotting the running time with varying number of corrupt worker nodes (up to <1/2 fraction).

2Storage redundancy is defined as the ratio of the size of the encoded matrix and the size of the raw data matrix.

Notation

We denote vectors by bold small letters (e.g., x, y, z, etc.)

And Matrices By Bold Capital Letters (E.G.,

A, F, S, X, etc.). We denote the amount of storage required by a matrix X by |X|.

For Any Positive

integer n ∈N, we denote the set {1, 2, . . , n} by [n]. For n1, n2 ∈N, where n1 ≤n2, we write [n1 : n2] to denote the set {n1, n1 + 1, . , n2}. For any vector u ∈Rn and any set S ⊂[n], we write uS to denote the |S|-length vector, which is the restriction of u to the coordinates in the set S. The support of a vector u ∈Rn is defined by supp(u) := {i ∈[n] : ui̸ = 0}. We say that a vector u ∈Rn is t-sparse if |supp(u)| ≤t.

While stating our results, we assume that performing the basic arithmetic operations (addition, subtraction, multiplication, and division) on real numbers takes unit time.

Problem Setting And Our Results

Given a dataset consisting of n labelled data points (xi, yi) ∈Rd × R, i ∈[n], we want to learn a model/parameter vector w ∈Rd, which is a minimizer of the following empirical risk minimization problem:

where fi(w), i = 1, 2, . . , n, denotes the risk associated with the i’th data point with respect to w and

Pn

i=1 fi(w) the average empirical risk associated with the n data points with respect to w. Our main focus in this paper is on generalized linear models (GLM), where fi(w) = ℓ(⟨xi, w⟩; yi) for some differentiable loss function ℓ. Here, each fi : Rd →R is differentiable, h : Rd →R is convex but not necessarily differentiable, and ⟨xi, w⟩is the dot product of xi and w. We do not necessarily need each fi to be convex, but we require f(w) to be a convex function. Note that f(w) + h(w) is a convex function. In the following we study different algorithms for solving (1) to learn a GLM.

Proximal Gradient Descent

We can solve (1) using Proximal Gradient Descent (PGD). This is an iterative algorithm, in which we choose an arbitrary/random initial w0 ∈Rd, and then update the parameter vector according to the following

where αt is the step size or the learning rate at the t’th iteration, determining the convergence behaviour. There are standard choices for it; see, for example, [BV04, Chapter 9]. For any h and α, the proximal

Observe that if h = 0, then proxh,α(w) = w for every w ∈Rd, and PGD reduces to the classical gradient descent (GD). This encompasses several important optimization problems related to learning, for which prox operator has a closed form expression; some of these problems are given below.

• Lasso. Here Fi(W) = 1

2(⟨xi, w⟩−yi)2 and h(w) = λ∥w∥1. It turns out that proxh,α(z) for Lasso is equal to the soft-thresholding operator Sλα(z) [Tib15], which, for j ∈[d], is defined as

Zj −Λα

if zj > λα. • SVM dual. Jaggi [Jag13] showed an equivalence between the dual formulation of Support Vector Machines (SVM) and Lasso. Hence, SVM dual is also a special case of (1).

• Constrained optimization. We want to solve a constrained minimization problem minw∈C f(w), where C ⊆Rd is a closed, convex set. Define an indicator function IC for C as follows: IC(w) := 0, if w ∈C; and IC(w) := ∞, otherwise. Now, observe the following equivalence

W∈C F(W) ⇐⇒Min

w∈Rd f(w) + IC(w). If we solve the RHS using PGD, then it can be easily verified that the corresponding proximal operator is equal to the projection operator onto the set C [Tib15]. So, the proximal gradient update step is to compute the usual gradient and then project it back onto the set C.

• Logistic regression. Here fi is the logistic function, defined as

,

where ui = ⟨xi, w⟩, and h = 0. As noted earlier, since h = 0, PGD reduces to GD for logistic regression.

Since Fi’S And H Are Differentiable,

we can alternatively solve this simply using GD. Let X ∈Rn×d denote the data matrix, whose i’th row is equal to the i’th data point xi. For simplicity,

M +J. In

a distributed setup, all the data is distributed among m worker machines (worker i has Xi) and master updates the parameter vector using the update rule (2). At the t’th iteration, master sends wt to all the workers; worker i computes the gradient (denoted by ∇if(wt)) on its local data and sends it to the master; master aggregates all the received m local gradients to obtain the global gradient

Now, master updates the parameter vector according to (2) and obtains wt+1. Repeat the process until convergence. If full gradients are too costly to compute. Updating the parameter vector in each iteration of PGD according to (2) requires computing full gradients. This may be prohibitive in large-scale applications, where each machine in a distributed framework has a lot of data, and computing full gradients at local machines may be too expensive and becomes the bottleneck. In such scenarios, there are two alternatives to reduce this per-iteration cost: (i) Coordinate Descent (CD), in which we pick a few coordinates (at random), compute the partial gradient along those, and descent along those coordinates only, and (ii) Stochastic Gradient Descent (SGD), in which we sample a data point at random, compute the gradient on that point, and descent along that direction. These are discussed in Section 2.2 and Section 6.1, respectively.

Oordinate Descent

For the clear exposition of ideas, we focus on the non-regularized empirical risk minimization from (1) (i.e., taking h = 0) for learning a generalized linear model (GLM). This can be generalized to objectives with (non-)differentiable regularizers [BKBG11, ST11]. Let X ∈Rn×d denote the data matrix and y ∈Rn the corresponding label vector. To make it distinct from the last section, we denote the objective function by φ and write it as φ(Xw; y) to emphasize that we want to learn a GLM, where the objective function depends on the data points only through their inner products with the parameter vector. Formally, we want to

For U ⊆[d], we write ∇Uφ(Xw; y) to denote the gradient of φ(Xw; y) with respect to wU, where wU denotes the |U|-length vector obtained by restricting w to the coordinates in U. To make the notation less cluttered, let φ′(Xw; y) denote the n-length vector, whose i’th entry is equal to ℓ′(⟨xi, w⟩; yi) :=

∂

∂uℓ(u; yi)|u=⟨xi,w⟩. Note that ∇φ(Xw; y) = XT φ′(Xw; y) and that ∇Uφ(Xw; y) = XT

Uφ′(Xw; Y), Where Xu Denotes The N×|U|

matrix obtained by restricting the column indices of X to the elements in U. Coordinate descent (CD) is an iterative algorithm, where, in each iteration, we choose a set of coordinates and update only those coordinates (while keeping the other coordinates fixed). In distributed CD, we take advantage of the parallel architecture to improve the running time of (centralized) CD. In the distributed setting, we divide the data matrix vertically into m parts and store the i’th part at the i’th worker node.

Concretely, assume, for simplicity, that m divides d. Let X = [X1 X2 . . Xm] and w = [wT

M Vector. Each Worker I Stores Xi And Is

responsible for updating (a few coordinates of) wi – hence the terminology, model-parallelism. We store the label vector y at the master node. In coordinate descent, since we update only a few coordinates in each round, there are a few options on how to update these coordinates in a distributed manner: Subset of workers: Master picks a subset S ⊂[m] of workers and asks them to update their wi’s [RT16].

This may not be good in the adversarial setting, because if only a small subset of workers are updating their parameters, the adversary can corrupt those workers and disrupt the computation.

Subset Of Coordinates For All Workers:

All the worker nodes update only a subset of the coordinates of their local parameter vector wi’s. Master can (deterministically or randomly) pick a subset U (which may or may not be different for all workers) of f ≤d/m coordinates and asks each worker to updates only those coordinates. If master picks U deterministically, it can cycle through and update all coordinates of the parameter vector in ⌈d/mf⌉iterations.

In Algorithm 1, we give the distributed CD algorithm with the second approach, where all worker nodes update the coordinates of their local parameter vectors for a single subset U. We will adopt this approach in our method to make the distributed CD Byzantine-resilient. Let r =

M. For Any I ∈[M], Let

wi = [wi1 wi2 . . wir]T and Xi = [Xi1 Xi2 . Xir], where Xij is the j’th column of Xi. For any i ∈[m] and U ⊆[r], let wiU denote the |U|-length vector that is obtained from wi by restricting its entries to the coordinates in U; similarly, let XiU denote the n × |U| matrix obtained by restricting the column indices of Xi to the elements in U.

In Algorithm 1, for each worker i to update wi according to (6), where the partial gradient of φ with

Iuφ′(Pm

j=1 Xjwj; y) and worker i has only (Xi, wi), every other worker j sends Xjwj to the master, who computes φ′(Pm

J=1 Xjwj; Y)5 And Sends It Back To All The

workers. Observe that, even if one worker is corrupt, it can send an adversarially chosen vector to make the computation at the master deviate arbitrarily from the desired computation, which may adversely affect the update at all the worker nodes subsequently.6 Similarly, corrupt workers can send adversarially chosen information to affect the stopping criterion.

3Here we are not optimizing the average of loss functions – since n is a fixed number, this does not affect the solution space. 4After the 1st iteration, worker i need not multiply Xi with wi to obtain Xiwi in every iteration; as only a few coordinates of wi are updated, it only needs to multiply those columns of Xi that corresponds to the updated coordinates of wi.

5Note that even after computing Xw, master needs access to the labels yi, i = 1, 2, . . , n to compute φ′(Xw; y). Since y ∈Rn is just a vector, we can either store that at master, or, alternatively, we can encode y distributedly at the workers and master can recover that using the method developed in Section 4 for Byzantine-resilient distributed matrix-vector multiplication, where the matrix is an identity matrix and vector is equal to y.

6Specifically, suppose the i’th worker is corrupt and the adversary wants master to compute φ′(Xw + e; y) for any arbitrary vector e ∈Rn of its choice, then the i’th worker can send Xiwi + e to the master.

Algorithm 1 Distributed Coordinate Descent

1: Initialize. Each worker i ∈[m] starts with an arbitrary/random wi ∈Rr, where r =

M And, For

simplicity, we assume that m divides d. 2: while (until the stopping criteria at master is not satisfied) do

:

Worker i computes Xiwi and sends it to the master node.4

:

Worker i receives (U ⊆[r], φ′(Xw; y)) from the master node.

:

Worker i updates its local parameter vector as (where ∇iUφ(Xw; y) = XT

while keeping the other coordinates of wi unchanged, and sends the updated wi to the master.

:

Master receives {Xiwi}i∈[m] from the m workers.

Aster First Computes Xw = Pm

i=1 Xiwi and then computes φ′(Xw; y).

:

Master picks U ⊆[r] (where U can be picked either randomly or in a round-robin fashion) and sends (U ⊆[r], φ′(Xw; y)) to all workers.

Adversary Model

We want to perform the distributed computation described in Section 2.1 and Section 2.2 under adversarial attacks, where the corrupt nodes may provide erroneous vectors to the master node. Our adversarial model is described next.

In our adversarial model, the adversary can corrupt at most t < m

Worker Nodes7, And The Compromised

then instead of sending the true vector, it may send an arbitrary vector to disrupt the computation. We refer to the corrupt nodes as erroneous or under the Byzantine attack. We can also handle asynchronous updates, by dropping the straggling nodes beyond a specified delay, and still compute the correct gradient due to encoding.

Therefore we treat updates from these nodes as being “erased”.

We Refer To These As

erasures/stragglers. For every worker i that sends a message to the master, we can assume, without loss of generality, that the master receives ui+ei, where ui is the true vector and ei is the error vector, where ei = 0 if the i’th node is honest, otherwise can be arbitrary. We assume that at most t nodes can be adversarially corrupt and at most s nodes can be stragglers, where s and t are some constants less than 1

That We Will

decide later. Note that the master node does not know which t worker nodes are corrupted (which makes this problem non-trivial to solve), but knows t. We propose a method that mitigates the effects of both of these anomalies.

Remark 1. A well-studied problem is that of asynchronous distributed optimization, where the workers can have different delays in updates [DB13]. One mechanism to deal with this is to wait for a subset of responses, before proceeding to the next iteration, treating the others as missing (or erasures) [KSDY17]. Byzantine attacks are quite distinct from such erasures, as the adversary can report wrong local gradients, requiring the master node to create mechanisms to overcome such attacks. If the master node simply aggregates the collected updates as in (4), the computed gradient could be arbitrarily far away from the true one, even with a single adversary [MGR18].

7Our results also apply to a slightly different adversarial model, where the adversary can adaptively choose which of the t worker nodes to attack at each iteration. However, in this model, the adversary cannot modify the local stored data of the attacked node, as otherwise, over time, it can corrupt all the data, making any defense impossible.

Our Approach To Gradient Computation

Recall that fi(w) = ℓ(⟨xi, w⟩; yi) for some differentiable loss function ℓ, and the gradient of fi at w is equal to ∇fi(w) = (xi)T ℓ′(⟨xi, w⟩; yi), where ℓ′(⟨xi, w⟩; yi) :=

∂Uℓ(U; Yi)|U=⟨Xi,W⟩. Note That ∇Fi(W) ∈Rd Is A

column vector. Let f ′(w) denote the n-length vector whose i’th entry is equal to ℓ′(⟨xi, w⟩; yi). With this

I=1 Fi(W), We Have ∇F(W) = 1

nXT f ′(w). Since n is a constant, it is enough to compute XT f ′(w). So, for simplicity, in the rest of the paper we write

A natural approach to computing the gradient ∇f(w) is to compute it in two rounds: (i) compute f ′(w) in the 1st round by first multiplying X with w and then master locally computes f ′(w) from Xw (master can do this locally, because Xw is an n-dimensional vector whose i’th entry is equal to ⟨xi, w⟩and (f ′(w))i = ℓ′(⟨xi, w⟩; yi));8 and then (ii) compute ∇f(w) = XT f ′(w) in the 2nd round by multiplying XT with f ′(w). So, the task of each gradient computation reduces to two matrix-vector (MV) multiplications, where the matrices are fixed and vectors may be different each time. To combat against the adversarial worker nodes, we do both of these MV multiplications using data encoding and real-error correction; see Figure 1 on page 17 for a pictorial description of our approach.

A two-round approach for gradient computation has been proposed for straggler mitigation in [LLP+18], but our method for MV multiplication differs from that fundamentally, as we have to provide adversarial protection. Note that in the case of stragglers/erasures we know who the straggling nodes are, but this infor- mation is not known in the case of adversarial nodes, and master needs to decode without this information in the context of Byzantine adversaries. This is slightly different from the standard error correcting codes (over finite fields) as the matrix entries in machine learning applications are from reals. In this case, we use ideas from real-error correction (or sparse reconstruction) from the compressive sensing literature [CT05], and using which we develop an efficient decoding at master, which also gives rise to our sparse encoding matrix; see Section 4 for more details. For decoding efficiently, we crucially leverage the block error pattern and design a decoding method at master, which, interestingly, requires just one application of the sparse recovery method on a vector of size m, the number of workers, which may be much smaller than the data dimensions n and d, thereby making the decoding computationally efficient. Our encoding matrix (given in (11), designed for MV multiplication) is very sparse and has a regular pattern of non-zero entries, which also makes it applicable for making coordinate-descent (CD) Byzantine-resilient. We emphasize that a general- purpose code for MV multiplication may not be applicable for CD, as each CD iteration requires updating only a few coordinates of the parameter vector, which makes it fundamentally different (and arguably more complicated to robustify) than GD iterations; see Section 3.2 and Section 5 for more details. Since iterative algorithms (such as GD and CD) require repeated parameter updates, it is crucial to have a method that has low computational complexity, both at the worker nodes as well as at the master node, and our coding solutions for both GD and CD achieve that, in addition to being highly storage efficient; see Theorem 1 for GD and Theorem 2 for CD.

Coming back to our two-round approach for gradient computations using MV multiplications, for the 1st round, we encode X using a sparse encoding matrix S(1) = [(S(1)

I X

at the i’th worker node; and for the 2nd round, we encode XT using another sparse encoding matrix

M )T ]T , And Store S(2)

i XT at the i’th worker node. Now, in the 1st round of the gradient computation at w, the master node broadcasts w and the i’th worker node replies with S(1)

I Xw

(a corrupt worker may report an arbitrary vector); upon receiving all the vectors, the master node applies error-correction procedure to recover Xw and then locally computes f ′(w) as described above. In the 2nd round, the master node broadcasts f ′(w) and similarly can recover XT f ′(w) (which is equal to the gradient) at the end of the 2nd round. So, it suffices to devise a method for multiplying a vector v to a fixed matrix A in a distributed and adversarial setting. Since this is a linear operation, we can apply error correcting codes over real numbers to perform this task. We describe it briefly below.

8Note that even after computing Xw, master needs access to the labels yi, i = 1, 2, . . , n to compute f′(w). See Footnote 5 for a discussion on how master can get access to the labels. A trivial approach. Take a generator matrix G of any real-error correcting linear code. Encode A as AT G =: B. Divide the columns of B into m groups as B = [B1 B2 . Bm], where worker i stores Bi.

Master broadcasts v and each worker i responds with vT Bi + eT

I , Where Ei = 0 If The I’Th Worker Is Honest,

otherwise can be arbitrary. Note that at most t of the ei’s can be non-zero. Responses from the workers can be combined as vT B + eT . Since every row of B is a codeword, vT B = vT AT G is also a codeword. Therefore, one can take any off-the-shelf decoding algorithm for the code whose generator matrix is G and obtain vT AT . For example, we can use the Reed-Solomon codes (over real numbers) for this purpose, which only incurs a constant storage overhead and tolerates optimal number of corruptions (up to < 1

). Note

that we need fast decoding, as it is performed in every iteration of the gradient computation by the master. As far as we know, any off-the-shelf decoding algorithm “over real numbers” requires at least a quadratic computational complexity, which leads to Ω(n2 + d2) decoding complexity per gradient computation, which could be impractical.

The trivial scheme does not exploit the block error pattern which we crucially exploit in our coding scheme to give a ∼O((n + d)m) time decoding per gradient computation, which could be a significant improvement our coding scheme enables a trade-off(in terms of storage and computation/communication overhead at the master and the worker nodes) with Byzantine adversary tolerance, without compromising the efficiency at the master node. We also want encoding to be efficient (otherwise it defeats the purpose of data encoding) and our sparse encoding matrix achieves that. Our main result for the Byzantine-resilient distributed gradient computation is as follows, which is proved in Section 4: Theorem 1 (Gradient Computation). Let X ∈Rn×d denote the data matrix. Let m denote the total number of worker nodes. We can compute the gradient exactly in a distributed manner in the presence of t corrupt worker nodes and s stragglers, with the following guarantees, where ϵ > 0 is a free parameter.

K

. • Total storage requirement is roughly 2(1 + ϵ)|X|. • Computational complexity for each gradient computation:

– At Each Worker Node Is O((1 + Ε) Nd

m ). – at the master node is O((1 + ϵ)(n + d)m). • Communication complexity for each gradient computation:

real numbers. – master broadcasts (n + d) real numbers.



. Remark 2. The statement of Theorem 1 allows for any s and t as long as (s + t) ≤

As We Are

handling both erasures and errors in the same way9 the corruption threshold does not have to handle s and t separately. To simplify the discussion, for the rest of the paper, we consider only Byzantine corruption, and denote the corrupted set by I ⊂[m] with |I| ≤t, with the understanding that this can also work with stragglers.

In Theorem 1, ϵ is a design choice and a free parameter that can take any value in the interval [0, m−1], where ϵ = 0 implies no corruption and ϵ = m −1 implies that corruption threshold t can be anything up to

M−1

2 . If we want to tolerate t corrupt workers, then ϵ must satisfy ϵ ≥

M−2T.10

9When there are only stragglers, one can design an encoding scheme where both the master and the worker nodes oper- ate oblivious to encoding, while solving a slightly altered optimization problem [KSDY17], in which gradients are computed approximately, leading to more efficient straggler-tolerant GD.

10We could have written everything in terms of t, m, n, d, but we chose to introduce another variable ϵ which, in our opinion, clearly brings out the tradeoffbetween the corruption threshold and the resource requirements without cluttering the expressions. Remark 3 (Comparison with the plain distributed PGD). We compare the resource requirements of our method with the plain distributed PGD (which provides no adversarial protection), where all the data points are evenly distributed among the m workers. In each iteration, master sends the parameter vector w to all the workers; upon receiving w, all workers compute the gradients on their local data in O( nd

M ) Time (Per

worker) and send them to the master; master aggregates them in O(md) time to obtain the global gradient and then updates the parameter vector using (2). In our scheme (i) the total storage requirement is O(1 + ϵ) factor more;11 (see also Remark 4) (ii) the amount of computation at each worker node is O(1 + ϵ) factor more; (iii) the amount of computation at the

Master Node Is O((1 + Ε)(1 + N

d )) factor more, which is comparable in cases where n is not much bigger than

D; (Iv) Master Broadcasts (1 + N

d ) factor more data, which is comparable if n is not much bigger than d; and



factor more data, which is O(1 + ϵ) – a constant factor – as long as n = O(dm). Remark 4. Let m be an even number. Note that we can get the corruption threshold t to be any number less than m/2, but at the expense of increased storage and computation. For any δ > 0, if we want to get δ close to m/2, i.e., t = m/2 −δ, then we must have (1 + ϵ) ≥m/2δ. In particular, at ϵ = 2, we can tolerate up to m/3 corrupt nodes, with constant overhead in the total storage as well as on the computational complexity.

Note that when δ is a constant, i.e., t is close to

If T = M−1

2 , then ϵ = m −1. In this case, our storage redundancy factor is O(m). In contrast, the trivial scheme (see “trivial approach” on page 9) does better in this regime and has only a constant storage overhead, but at the expense of an increased decoding complexity at the master, which is at least quadratic in the problem dimensions d and n, whereas, our decoding complexity at the master always scales linearly with d and n. If we always want a constant storage redundancy for all values of the corruption threshold t, we can use our

Coding Scheme If T ≤C · M−1

2 , where c < 1 is a constant, and use the trivial scheme if t is close to m−1 2 .

Time. Note That O(Nd) Is Equal To The Time

required for distributing the data matrix X among m workers (for running the distributed gradient descent algorithms without the adversary); and the encoding time in our scheme (which results in an encoded matrix that provides Byzantine-resiliency) is a factor of (2t + 1) more.

Remark 5. Our scheme is not only efficient (both in terms of computational complexity and storage re-

Quirement), But It Can Also Tolerate Up To ⌊M−1

2 ⌋corrupt worker nodes (by taking ϵ = m −1 in Theorem 1). It is not hard to prove that this bound is information-theoretically optimal, i.e., no algorithm can tolerate

⌈M

2 ⌉corrupt worker nodes, and at the same time correctly computes the gradient.

Our Approach To Coordinate Descent

We use data encoding and add redundancy to enlarge the parameter space. Specifically, we encode the data matrix X using an encoding matrix R = [R1 R2 . . Rm], where each Ri is a d×p matrix (with pm ≥d), and store XRi at the i’th worker. Define eXR := XR. Now, instead of solving (5), we solve the encoded problem arg minv∈Rpm φ( eXRv; y) using Algorithm 1 (together with decoding at the master); see Figure 2 on page 25 for a pictorial description of our algorithm. We design the encoding matrix R such that at every iteration of our algorithm, updating any (small) subset of coordinates of vi’s (let v = [vT

M]) Automatically

updates some (small) subset of coordinates of w; and, furthermore, by updating those coordinates of vi’s, we can efficiently recover the correspondingly updated coordinates of w, despite the errors injected by the adversary. In fact, at any iteration t, the encoded parameter vector vt and the original parameter vector wt satisfies vt = R+wt, where R+ := RT (RRT )−1 is the Moore-Penrose pseudo-inverse of R, and wt evolves in the same way as if we are running Algorithm 1 on the original problem.

11For example, by taking ϵ = 2, our method can tolerate m/3 corrupt worker nodes. So, we can tolerate linear corruption with a constant overhead in the resource requirement, compared to the plain distributed gradient computation which does not provide any adversarial protection.

We will be effectively updating the coordinates of the parameter vector w in chunks of size (m −2t) or its integer multiples (where t is the number of corrupt workers). In particular, if each worker i updates k coordinates of vi, then k(m −2t) coordinates of w will get updated. For comparison, Algorithm 1 updates km coordinates of the parameter vector w in each iteration, if each worker updates k coordinates in that iteration.

As described in Algorithm 1 for the Byzantine-free CD, in order to update its local parameter vector wi according to (6), worker i needs access to φ′(Xw; y), which master computes after receiving {Xjwj}j∈[m] from the workers. In our Byzantine-resilient algorithm for CD also master will need to compute Xw in every CD iteration, and for this purpose, we employ the same encoding-decoding procedure for MV multiplication that we used in the first round of gradient computation, as described in Section 2.4. In particular, to make the notation distinct from gradient computation, in order to compute Xw, we encode X using an encoding

T

m]T , where each Li is a p′×n matrix (with p′m ≥n) and worker i stores eXL i = LiX. Note that in order to compute Xw, in the first round of gradient computation as described in Section 2.4, master broadcasts w to all the workers and each worker i computes eXL

I W And Sends It The The Master (Corrupt

workers may report arbitrary vectors), who then decodes and obtains Xw. However, in coordinate descent, though master wants to compute Xw in each CD iteration, we can significantly improve the computation required at each worker: since only a few coordinates of the original parameter vector w are updated in each CD iteration, master needs to send only those updated coordinates, and workers need to preform MV multiplication with a much smaller matrix, whose number of columns is equal to the number of updated coordinates of w that they receive from master. Thus, the computational complexity in each CD iteration at worker is proportional to the number of coordinates updated in each CD iteration, as desired.

Our main result for the Byzantine-resilient distributed coordinate descent is stated below, which is proved in Section 5. Theorem 2 (Coordinate Descent). Under the setting of Theorem 1, our Byzantine-resilient distributed CD algorithm has the following guarantees, where ϵ > 0 is a free parameter.

K

. • Total storage requirement is roughly 2(1 + ϵ)|X|. • If each worker i updates τ coordinates of vi, then

Τm

1+ϵ coordinates of the corresponding w gets updated.

– The Computational Complexity In Each Iteration

∗at each worker node is O(nτ). ∗at the master node is O((1 + ϵ)nm + τm2).



. Remark 6 (Comparison with the plain distributed CD). We compare the resource requirements of our method with the plain distributed CD described in Algorithm 1 that does not provide any adversarial protec- tion. Let ϵ be any number in the interval [0, m −1] – for illustration, we can take ϵ = 2, which means t ≤m workers are corrupt. In Algorithm 1, if each worker i updates

+Ε Coordinates

of w) in each iteration, then (i) each worker requires O( nτ

Wi; (Ii) Master Requires O(Nm) Time To Compute Pm

i=1 Xiwi from {Xiwi}i∈[m]; (iii) each worker sends n real numbers (required for Xiwi) to master; and (iv) master broadcasts n real numbers (required for φ′(Xw; y)). In our scheme (i) the total storage requirement is O(1+ϵ) factor more; (ii) the amount of computation at each worker node is O(1+ϵ) factor more; (iii) the amount of computation at the master node is O((1+ϵ)+ τm

N )

factor more – typically, since τ is a constant and number of workers is much less than n, this again could be



factor more data, which could be a constant if τm is smaller



factor more data, where the 1st term is much smaller than 1 as τ is typically a constant, and the 2nd term is close to zero as (1 + ϵ) is always upper-bounded by m.

Remark 7 (Comparison with the replication-based strategy). One simple way to make Algorithm 1 Byzantine- resilient is using repetition code, where we first divide the set of m workers into

T+1 Groups Of Size (2T + 1)

each and also divide the data matrix as X = [X1 X2 . . X

M

2t+1 ] (assume, for simplicity, that (2t+1) divides m). Now, store the i’th block Xi at the (2t + 1) workers in the i’th group of workers. Let the parameter

Wtm

2t+1 ]T . In each CD iteration, the local parameter updates in any wi is replicated at (2t + 1) different workers in the i’th group of workers, and since at most t workers are corrupt, master can do a majority vote for decoding. Note that the total storage and the computation at workers in this scheme grow linearly by a factor of (2t+1), where t is the number of corruption, which could be significant. In contrast, the method that we propose can tolerate linear corruption, say, t = m

, With A

constant overhead in storage and computational complexity. The Remarks 2, 4, 5 are also applicable for Theorem 2.

Related Work

There has been a significant recent interest in using coding-theoretic techniques to mitigate the well-known straggler problem [DB13], including gradient coding [TLDK17,RTDT18,CP18,HRSH18], encoding compu- tation [LLP+18, DCG16, DCG19], and data encoding [KSDY17, KSDY19]. However, one cannot directly apply the methods for straggler mitigation to the Byzantine attacks case, as we do not know which up- dates are under attack. Distributed computing with Byzantine adversaries is a richly investigated topic since [LSP82], and has received recent attention in the context of large-scale distributed optimization and learning [BMGS17, CSX17, CWCP18, YCRB18, AAL18, SX19, XKG19, YCRB19, GV19, RWCP19, LXC+19, GHYR19, YLR+19, DD20b, DD20a, HKJ20].

These can be divided into three categories: (i) One which assume explicit statistical models for data across workers (e.g., data drawn i.i.d. from a probability distri- bution) and analyze gradient descent [CSX17, YCRB18, SX19, YCRB19, GHYR19]. (ii) Other set of works make no probabilistic assumption on data, and optimize through stochastic methods (e.g., stochastic gra- dient descent) [BMGS17, AAL18, GV19, XKG19, LXC+19, RWCP19, DD20a, DD20b, HKJ20] and also with deterministic methods (e.g., gradient descent) [DD20a, DD20b]. Note that none of these two sets of works do data encoding and work with data as it is, and provide Byzantine resilience by applying some robust aggregation procedures (e.g., geometric median, coordinate-wise median, outlier-filtering, etc.) at the mas- ter for aggregating gradients. (iii) Another line of work which is most relevant to ours provide Byzantine resiliency using redundant computations, either by encoding the gradients [CWCP18] or by encoding the data itself [YLR+19]. Note that [RWCP19] combines both redundant computations and do a hierarchical robust aggregation and not is directly comparable to ours.

Note that the statistical nature of data/analysis in the first two sets of works leads to a statistical approximation error in the convergence rates, which is also intensified by the inaccuracy of the robust gradient aggregation procedure. One of the main focuses in these works is typically on obtaining faster convergence (where the goal is to match the convergence rate of plain SGD/GD) and as good an approximation error as possible. Note that the approximation error in all these works scales at least as Ω(

D), Where D Is The

dimension of the model parameter vector, which may be significant in high-dimensional settings. Moreover, need to make some assumptions on the data, and furthermore, master has to apply a non-trivial decoding for gradient aggregation, which requires significantly more time than what our decoding requires.

For

example, filtering-based decoding [SX19, DD20a, DD20b], median-based decoding [CSX17, YCRB18], and heuristic approaches [BMGS17], all have a super-linear complexity in m – in fact, the filtering-based method as in [SX19, DD20a, DD20b] (which is the most effective in terms of the approximation error) requires O(m3d) time. In contrast, our decoding has a linear dependence on both m and d. Note that, unlike the first two categories, the third line of work (to which ours also belongs) gives deterministic guarantees and work with arbitrary datasets, with no probabilistic assumptions; we elaborate on these and do a detailed comparison with ours below. We skip the comparison with the first two categories, as it would not be a fair comparison because the underlying setting is different – results in the first two categories are based on statistical assumptions on data/algorithm and inaccurate gradient recovery, whereas, results in the third category make no assumption on the data/algorithm and allow exact gradient recovery.

We want to emphasize that all these works use gradient descent (GD) or stochastic gradient descent (SGD) as their optimization algorithm, which is a data-parallelization method; in this paper, additionally, we also use coordinate descent (CD) algorithm for optimization, which is a model-parallelization method and is preferred over GD in some applications; see Section 1 for more details on this. As will be evident from Section 5, making CD secure against Byzantine attacks is arguably more intricate than securing GD.

We divide this section into three categories: first we compare the redundancy-based methods for GD in Section 3.1, and then CD in Section 3.2. Since we use matrix-vector (MV) multiplication as a core subroutine for both GD and CD, we also compare related work on this in Section 3.3.

Gradient Descent (Gd)

In this section, we do a detailed comparison with [CWCP18] and [YLR+19], which are the closest related works that also combat Byzantine adversaries using redundant computations.

Workers Are Corrupt. The Coding Scheme Of Chen Et

al. [CWCP18], which they called Draco, requires repetition of each data point (2t + 1) times, storing each copy at different workers. This gives the storage redundancy factor of (2t + 1) in Draco, whereas, our coding method requires storage redundancy factor of 2(1 + ϵ) =

Constant (< 1

each GD iteration (than simply computing the gradients as in plain distributed GD), the computational cost at workers also grows by the same factor, which is a significant downside of their scheme. In contrast, our

M

m−2t) more computation at worker, which is a constant even if t is a constant (< 1

)

fraction of m. This significantly reduces the computation time at the worker nodes in our scheme compared to Draco, without sacrificing much on the computation time required by the master node – the decoding at master in Draco takes O(md) time, whereas, our scheme requires O(

M−2T(1 + N

d )) more than Draco. In high-dimensional settings, where n is not much bigger than d, and

T Is A Constant (< 1

2) fraction of m, this overhead is constant. Overall, for a constant fraction of corruption,

Say, T = M

3 , Draco requires Ω(t) times more storage and computation at workers than our scheme (which could be significant in large-scale settings), and requires Ω(1 + n

D ) Times Less Computation At Master. Note

that the computation time at workers scales at least as Ω( nd

M ), Which Dominates The Time Taken By Master

(since n, d are typically much larger than m), so our scheme will be faster than Draco with respect to the overall running time. Note that the coding in Draco is restricted to data replication redundancy, as they encode the gradient as done in [TLDK17], enabling application to (non)-convex problems; in contrast, we encode the data enabling significantly smaller redundancy, and apply it to learn generalized linear models, and is also applicable to MV multiplication.

12To highlight the storage redundancy gain of our method over that of Draco, consider the following two concrete scenarios, where the data matrix X ∈Rn×d consists of nd real numbers: (i) In a large setup with m = 1000 worker nodes, if we want resiliency against t = 100 corrupt nodes (1/10 nodes are corrupt), our method requires redundancy of 2.5, whereas Draco requires redundancy of 201 (i.e., we need to store only 2.5 × nd real numbers, whereas Draco stores 201 × nd real numbers), a multiplicative-factor of > 80 more than ours. (ii) In a moderate setup with m = 150 and t = 50 (1/3 nodes are corrupt), the redundancy of our method is 6, whereas Draco requires redundancy of 101, a multiplicative-factor of ≈17 more than ours.

Yu et al. [YLR+19] (which is a concurrent work13) proposes Lagrange coded computing in a distributed framework to compute any multivariate polynomial of the input data and simultaneously provides resilience against stragglers, security against adversaries, and privacy of the dataset against collusion of workers. They leverage the Lagrange polynomial to create computation redundancy among workers, and using standard Reed-Solomon decoding, they can tolerate both erasures/stragglers and errors/adversaries. Their method provide privacy by adding random elements from the field (which in the case of gradient computation is the method in Shamir secret sharing scheme [Sha79] that is widely used in information-theoretically secure MPC protocols [CDN15] to provide privacy of users’ data. For the sake of comparison of the resource requirements of our scheme and the one in [YLR+19], consider the task of linear regression (the concrete machine learning application studied in [YLR+19]). In the following, we assume that m−1

M

1+2δ −1 in our setting; here δ can take any value in [0 : m−1

M

δ+1, which is roughly the same as ours. For example,

Corrupt Workers (I.E., Δ = M−3

6 ), the storage overhead of our scheme and of [YLR+19] is a

Multiplicative Factor Of 6 And

1+3/m ≈6, respectively. (ii) The encoding time complexity of our scheme is

O(Nd(M−2Δ)), Whereas, It Is O(M Log2(M) Nd

δ+1) in [YLR+19]. Note that for constant δ (i.e., corruption close to 1/2), the encoding time of our scheme is much less (by a factor of O(m log2(m))) than that of [YLR+19],

Log2(M))-Factor Less Time In

encoding than ours. (iii) The computation time at each worker per gradient computation in both our scheme and [YLR+19] is roughly the same – ours requires O(

+Δ) Time. (Iv) The

decoding time complexity per gradient computation in [YLR+19] is O(m log2(m)d), whereas, ours requires O((1 + ϵ)(n + d)m) time. Note that when n is not much bigger than d and we want a constant fraction of

Corruption, Say, M

3 corruption, then their decoding complexity is worse than ours by a logarithmic factor. Also note that our decoding algorithm is arguably simpler than theirs. (v) For per gradient computation,

N+D

1+2δ and d real numbers in ours and the scheme in [YLR+19]. Note that if n ≤dm and to tolerate a constant fraction of corruption, say, m

Corruption, Each Worker Sends Roughly

O(m) less data in our scheme than that of [YLR+19]. Overall, if we want tolerance against m

Corrupt Worker

nodes, then both our scheme and the one in [YLR+19] have similar resource requirements, except for that our scheme has a much better communication complexity (by a factor of O(m)) from workers to the master, whereas, the encoding time complexity (which is a one-time process) of [YLR+19] is better than ours by a

Oordinate Descent (Cd)

Even for the straggler problem, we are only aware of one work by Karakus et al. [KSDY19] that, in addition to distributed GD, also studies distributed CD, and that for quadratic problems (e.g., linear/ridge regression) only. It also does data encoding and achieves low redundancy and low complexity, by allowing convergence to an approximate rather than exact solution. As far as we know, ours is the first work that studies distributed CD under Byzantine attacks and provides an efficient solution, much better than the replication-based solution (see Remark 7).

At the heart of our solution for CD is the matrix-vector (MV) multiplication procedure that we develop in this paper; and it is the specific regular structure of our encoding matrix (given in (11), designed for the MV multiplication) that allows for partially updating the coordinates of the parameter vector in each CD iteration. Note that a general-purpose encoding matrix for MV multiplication may not be applicable for the CD algorithm.

It has been observed earlier in several works (see, for example, [LLP+18,DCG16]) that gradient compu- tation in GD for linear regression can be reduced to MV multiplication, and any general-purpose code for MV multiplication can be used to provide a solution for gradient computation. As far as we know, ours is the first paper that makes the connection of CD and MV multiplication, and provides an efficient solution 13Yu et al. [YLR+19] is concurrent to our conference versions in Allerton 2018 [DSD18] and ISIT 2019 [DSD19, DD19], on which this paper is based.

for CD (which is also resilient to Byzantine attacks) for learning generalized linear models. Note that, unlike GD, not any general-purpose code for MV multiplication can be used for CD: the main challenge in CD comes from the fact that we only update a small number of coordinates of the parameter vector in each CD iteration; when we encode the data and iteratively update some coordinates of the (encoded) parameter vector using the encoded data, we need to make sure that this update in the encoded parameter vector is reconciled with the update in the original parameter vector. This is fundamentally different from GD iterations. See Section 5 for more details.

Atrix-Vector Multiplication

For the task of a more fundamental problem of matrix-vector (MV) multiplication in the presence of Byzan- tine adversaries, which is at the core of the optimization algorithms in this paper, we are only aware of two concurrent works [YLR+19] (see Footnote 13) and [DCG19]14 that provide (coding-theoretic) solutions to this problem. In the following, we do a detailed comparison of our solution with both of these works and also discuss the (dis)similarities.

We have already done a detailed comparison with Yu et al. [YLR+19] (concurrent work, see Footnote 13) with respect to gradient descent in Section 3.1. For the problem of MV multiplication, the storage require- ment, computation time per worker, and communication complexity to/from workers is the same in both ours and [YLR+19]. The comparison of encoding time complexity is same as above; however, for a constant

Corruption, Say, M

3 corrupt workers, our method outperforms the one in [YLR+19] in terms of the decoding time complexity by a factor of O(log2(m)). Note that, unlike [YLR+19], we make a fundamental connection of handling Byzantine errors with the sparse reconstruction (or the real-error correction) problem from the compressive sensing literature [CT05].

Dutta et al. [DCG19] (concurrent work, see Footnote 14) focuses on matrix-vector (MV) multiplication. Though their main focus is on providing resilience against stragglers, they also mention that handling stragglers is very different than handling errors, as it requires to correct errors over real numbers, and, unlike stragglers, we do not know which workers are corrupt. Similar to our observation, they also note that since the matrices and vectors have entries from real numbers, the decoding problem reduces to the sparse reconstruction problem from the compressive sensing literature [CT05] and they also provide such a reduction. Apart from these similarities, our solution for MV multiplication differs from that of [DCG19] in several important ways: (i) [DCG19] provides a detailed solution to the distributed MV multiplication for the straggler problem for the case when the number of rows in the matrix is smaller than the number of workers nodes. As mentioned in [DCG19], this method can be easily generalized to the more general case when the matrix is of arbitrary dimension, in which case, first we can divide the rows of the matrix into several sub-matrices, each having number of rows smaller than the number of workers, and then apply the above method independently to each sub-matrix. This simple extension may work (without losing efficiency) for the straggler/erasure problem, however, leads to a highly inefficient solution for the adversary/error problem.

The reason being that, in the presence of Byzantine workers, if we solve the sparse reconstruction problem for each sub-matrix separately, this would be inefficient, as the decoding would then be computationally expensive. To remedy this, we exploit the block error pattern and use a simple idea of linearly combining the response vectors from each worker using coefficients drawn from an absolutely continuous distribution, so that we only need to do just one computation for solving the sparse construction problem. This significantly reduces the decoding complexity; see Section 4.1 for details.

(Ii) [Dcg19] Only Shows A Connection To

the sparse recovery problem, whereas, we provide a complete solution, with a concrete sparse recovery (or real-error correction) matrix and resource (encoding/decoding time, storage, communication) requirement analysis. (iii) Our encoding matrix (given in (11)) to encode data matrices of arbitrary dimensions is very sparse and highly structured which allows us to apply that construction to CD algorithm, which, as far we know, has not been connected with MV multiplication before. Also, ours is the first paper that provides a non-trivial and efficient (data encoding) solution to CD in the presence of a Byzantine adversary. (iv) We 14The conference version [DCG16] only studies the straggler problem, and the journal version [DCG19] briefly mentions how their results from [DCG16] can be extended to handle adversarial nodes, and we describe that in this section.

also want to mention that the focus in [DCG19] is on making the encoded matrix sparse (at the expense of increased computation at workers) so that workers need to compute shorter dot products, whereas, in this paper, we make the encoding matrix sparse (much sparser than the encoded matrix of [DCG19]) to get efficient encoding/decoding.

Our Solution To Gradient Computation

In this section, we describe the core technical part of our two-round approach for gradient computation described in Section 2.4 – a method for performing matrix-vector (MV) multiplication in a distributed manner in the presence of a malicious adversary who can corrupt at most t of the m worker nodes. Here, the matrix is fixed and we want to right-multiply a vector with this matrix.

Given a fixed matrix A ∈Rnr×nc and a vector v ∈Rnc, we want to compute Av in a distributed manner in the presence of at most t corrupt worker nodes; see Section 2.3 for details on our adversary model. Our method is based on data encoding and error correction over real numbers, where the matrix A is encoded and distributed among all the worker nodes, and the master node recovers the MV product Av using real-error correction; see Figure 1. We will think of our encoding matrix as S = [ST

M], Where Each Si Is

a p × nr matrix and pm ≥nr. We will derive the matrix S in Section 4.2. For the value of p, looking

, We Would Have P = 3N

m ). For i ∈[m], we store the matrix SiA at the i’th worker node. As described in Section 2, the computation proceeds as follows: The master sends v to all the worker nodes and

Receives {Siav + Ei}M

i=1 back from them. Let ei = [ei1, ei2, . . , eip]T for every i ∈[p]. Note that ei = 0 if the i’th node is honest, otherwise can be arbitrary. In order to find the set of corrupt worker nodes, master

Equivalently Writes {Siav + Ei}M

i=1 as p systems of linear equations.

where, for every i ∈[p], ˜ei = [e1i, e2i, . . , emi]T , and ˜Si is an m × nr matrix whose j’th row is equal to the i’th row of Sj, for every j ∈[m]. Note that at most t entries in each ˜ei are non-zero. Observe that

I=1 And {˜Siav + ˜Ei}P

i=1 are equivalent systems of linear equations, and we can get one from the other. Note that ˜Si’s constitute the encoding matrix S, which we have to design. In the following, we will design these matrices ˜Si’s (which in turn will determine the encoding matrix S), with the help of another matrix F, which will be used to find the error locations, i.e., identities of the compromised worker nodes. We will design the matrix F (of dimension k × m, where k < m – here k is determined by the error-correction capability, and we will set k = 2t; see Section 4.4 for more details) and the matrices ˜Si’s such that C.1 F˜Si = 0 for every i ∈[p].

C.2 For any t-sparse u ∈Rm, we can efficiently find all the non-zero locations of u from Fu. C.3 For any T ⊂[m] such that |T | ≥(m −t), let ST denote the |T |p × nr matrix obtained from S by restricting it to all the Si’s for which i ∈T . We want ST to be of full column rank.

If we can find such matrices, then we can recover the desired MV multiplication Av exactly: briefly, C.1 and C.2 will allow us to locate the corrupt worker nodes; once we have found them, we can discard all the information that the master node had received from them. This will yield ST Av, where ST is the |T |p × nr matrix obtained from S by restricting it to Si’s for all i ∈T , where T is the set of all honest worker nodes.

Now, by C.3, since ST is of full column rank, we can recover Av from ST Av exactly. Details follow. Suppose we have matrices F and ˜Si’s such that C.1 holds. Now, multiplying (8) by F yields

for every i ∈[p], where ∥˜ei∥0 ≤t. In Section 4.1, we give our approach for finding all the corrupt worker nodes with the help of any error locator matrix F. Then, in Section 4.2, we give a generic construction for

W ←−Proxh,Α(W −Α∇F(W))

Figure 1 This figure shows our 2-round approach to the Byzantine-resilient distributed gradient descent to optimize (1) for learning a generalized linear model. Since the gradient at w is equal to ∇f(w) = XT f′(w) (see (7)), we compute it in 2 rounds, using a matrix-vector (MV) multiplication as a subroutine in each round. In the 1st round, first we compute Xw, and then compute f′(w) from Xw – since the j’th entry of Xw is equal to ⟨xj, w⟩, we can compute f′(w) from Xw (see Section 2.4).

In the 2nd round we compute XT f′(w) – which is equal to ∇f(w) – using another application of MV multiplication. For a matrix A and a vector v, to make our distributed MV multiplication Av Byzantine-resilient, we encode A using a sparse matrix

M . . . St

m]T and distribute SiA to worker i (denoted by Wi). Note that in the first round, we have A = X, v = w, and we encode X using S(1), and in the second round, we have A = XT , v = f′(w), and encode XT using S(2). The adversary can corrupt at most t workers (the compromised ones are denoted in red color), potentially different sets of t workers in different rounds. The master node (denoted by M) broadcasts v to all the workers. Each worker performs the local MV product and sends it back to M. If Wi is corrupt, then it can send an arbitrary vector. Once the master has received all the vectors (out of which t may be erroneous), it sends them to the decoder (denoted by Dec), which outputs the correct MV product Av.

designing ˜Si’s (and, in turn, our encoding matrix S) such that C.1 and C.3 hold. In Section 4.3, we show how to compute the desired matrix-vector product Av efficiently, once we have discarded all the data from the corrupt works nodes. Then, in Section 4.4, we will give details of the error locator matrix F that we use in our construction.

Remark 8. As we will see in Section 4.2, the structure of our encoding matrix S is independent of our error locator matrix F. Specifically, the repetitive structure of the non-zero entries of S as well as their locations will not change irrespective of what the F matrix is. This makes our construction very generic, as we can choose whichever F suits our needs the best (in terms of how many erroneous indices it can locate and with what decoding complexity), and it won’t affect the structure of our encoding matrix at all – only the non-zero entries might change, neither their repetitive format, nor their locations!

Finding The Corrupt Worker Nodes

Observe that supp(˜ei) may not be the same for all i ∈[p], but we know, for sure, that the non-zero locations in all these error vectors occur within the same set of t locations. Let I = Sp

I=1 Supp(˜Ei), Which Is The Set

of all corrupt worker nodes. Note that |I| ≤t. We want to find this set I efficiently, and for that we note the following crucial observation. Since the non-zero entries of all the error vectors ˜ei’s occur in the same set I, a random linear combination of ˜ei’s has support equal to I with probability one, if the coefficients of the linear combination are chosen from an absolutely continuous probability distribution. This idea has appeared before in [ME08] in the context of compressed sensing for recovering arbitrary sets of jointly sparse signals that have been measured by the same measurement matrix.

Definition 1. A probability distribution is called absolutely continuous, if every event of measure zero occurs with probability zero. It is well-known that a distribution is absolutely continuous if and only if it can be represented as an integral over an integrable density function [Bil95, Theorem 31.8, Chapter 6]. Since Gaussian and uniform distributions have an explicit integrable density function, both are absolutely continuous. Conversely, dis- crete distributions are not absolutely continuous. Now we state a lemma from [ME08] that shows that a random linear combination of the error vectors (where coefficients are chosen from an absolutely continuous distribution) preserves the support with probability one.

I=1 Αi˜Ei, Where Αi’S Are Sampled I.I.D. From An

absolutely continuous distribution. Then with probability 1, we have supp(ˆe) = I. From (9) we have fi = F˜ei for every i ∈[p]. Take a random linear combination of fi’s with coefficients αi’s chosen i.i.d. from an absolutely continuous distribution, for example, the Gaussian distribution. Let

I=1 Αi˜Ei. Note That, With Probability

1, supp(˜e) is equal to the set of all corrupt worker nodes, and we want to find this set efficiently. In other words, given F˜e, we want to find supp(˜e) efficiently. For this, we need to design a k × m matrix F (where k < m) such that for any sparse error vector e ∈Rm, we can efficiently find supp(e) from f = Fe. Many such matrices have been known in the literature that can handle different levels of sparsity with varying decoding complexity. We can choose any of these matrices depending on our need, and this will not affect the design of our encoding matrix S. In particular, we will use a k × m Vandermonde matrix along with the Reed-Solomon type decoding, which can correct up to k/2 errors and has decoding complexity of O(m2); see Section 4.4 for details.

Time required in finding the corrupt worker nodes.

The Time Taken In Finding The Corrupt Worker

nodes is equal to the sum of the time taken in the following 3 tasks. (i) Computing F˜ei for every i ∈[p]: Note that we can get F˜ei by multiplying (8) with F. Since F is a k × m matrix, and we compute F˜hi(v) for p systems, this requires O(pkm) time. (ii) Taking a random linear combination of p vectors each of length m, which takes O(pm) time. (iii) Applying Lemma 2 (in Section 4.4) once to find the error locations, which takes O(m2) time. Since p is much bigger than m, the total time complexity is O(pkm).

Esigning The Encoding Matrix S

Now we give a generic construction for designing ˜Si’s such that C.1 and C.3 hold. Fix any k × m matrix F such that we can efficiently find e from Fe, provided e is sufficiently sparse. We can assume, without loss of generality, that F has full row-rank; otherwise, there will be redundant observations in Fe that we can discard and make F smaller by discarding the redundant rows. Let N(F) ⊂Rm denote the null-space of F. Since rank(F) = k, dimension of N(F) is q = (m −k). Let {b1, b2, . . , bq} be a basis of N(F), and let bi = [bi1 bi2 . bim]T , for every i ∈[q]. We set bi’s the columns of the following matrix F⊥:

The following property of F⊥will be used for recovering the MV product in Section 4.3. Claim 1. For any subset T ⊂[m], such that |T | ≥(m −t), let F⊥

The Restriction Of F⊥To The Rows In T . Then F⊥

T is of full column rank. Proof. Note that q = m −k, where k = 2t. So, if we show that any q rows of F⊥are linearly independent, then, this in turn will imply that for every T ⊂[m] with |T | ≥(m −t), the sub-matrix F⊥

T Will Have Full

column rank. In the following we show that any q rows of F⊥are linearly independent. To the contrary, suppose not; and let T ′ ⊂[m] with |T ′| = q be such that the q ×q matrix F⊥

T ′ Is Not A Full Rank Matrix. This

implies that there exists a non-zero c′ ∈Rq such that F⊥ T ′c′ = 0. Let b = F⊥c′. Note that b̸ = 0 (because columns of F⊥are linearly independent) and also that ∥b∥0 ≤m −q = k. Now, since FF⊥= 0, we have Fb = 0, which contradicts the fact that any k columns of F are linearly independent.

Now we design ˜Si’s. For i ∈[p], we set ˜Si as follows:



where l = q if i < p; otherwise l = nr −(p −1)q. The first (i −1)q and the last nr −[(i −1)q + l] columns of ˜Si are zero. This also implies that the number of rows in each Si is p = ⌈nr/q⌉. Claim 2. For every i ∈[p], we have F˜Si = 0.

Proof. By construction, the null-space of F is N(F) = span{b1, b2, . . , bq}, which implies that Fbi = 0, for every i ∈[q]. Since all the columns of ˜Si’s are either 0 or bj for some j ∈[q], the claim follows. The above constructed matrices ˜Si’s give the following encoding matrix Si for the i’th worker node:

All the unspecified entries of Si are zero. The matrix Si is for encoding the data for worker i. By stacking up the Si’s on top of each other gives us our desired encoding matrix S. To get efficient encoding, we want S to be as sparse as possible. Since S is completely determined by F⊥, whose columns are the basis vectors of N(F), it suffices to find a sparse basis for N(F). It is known that finding the sparsest basis for the null-space of a matrix is NP-hard [CP86]. Note that we can always find the basis vectors of N(F) by reducing F to its row-reduced-echelon-form (RREF) using the Gaussian elimination [HK71]. This will result in F⊥whose last q rows forms a q × q identity matrix. Note that q = m −k, where k = 2t. So, if the corruption threshold t is very small as compared to m, the F⊥that we obtain by the RREF will be very sparse – only the first 2t rows may be dense. Since computing S is equivalent to computing F⊥, and we can compute F⊥in O(k2m) time using the Gaussian elimination, the time complexity of computing S is also O(k2m).

Now we prove an important property of the encoding matrix S that will be crucial for recovery of the desired matrix-vector product. Claim 3. For any T ⊂[m] such that |T | ≥(m −t), let ST denote the |T |p × nr matrix obtained from S by restricting it to all the blocks Si’s for which i ∈T . Then ST is of full column rank.

Proof. For i ∈[p −1], let Bi = [(i −1)q + 1 : iq] and Bp = [(p −1)q + 1 : nr −(p −1)q], where we see Bi’s as a collection of some column indices. Consider any two distinct i, j ∈[p]. It is clear that for any two vectors u1 ∈Bi, u2 ∈Bj, we have supp(u1) ∩supp(u2) = φ, which means that all the columns in distinct Bi’s are linearly independent. So, to prove the claim, we only need to show that the columns within the same Bi’s are linearly independent. Fix any i ∈[p], and consider the |T |p × q sub-matrix S(i)

Authors:

Peder EZ Larson 1, 2,* , Jenna ML Bernard1, James A Bankson 3, Nikolaj Bøgh 4, Robert A Bok1, Albert P. Chen 5, Charles H Cunningham 6,7, Jeremy Gordon1, Jan-Bernd Hövener 8, Christoffer Laustsen 4, Dirk Mayer 9,10, Mary A McLean11 12, Franz Schilling13, James Slater1, Jean-Luc Vanderheyden5, 14, Cornelius von Morze 15, Daniel B Vigneron1, 2, Duan Xu1, 2, and the HP 13C

94143, Usa.

Denmark. 5 GE Healthcare, Menlo Park, California, USA. 6 Physical Sciences, Sunnybrook Research Institute, Toronto, Ontario, Canada.

8 Section Biomedical Imaging, Molecular Imaging North Competence Center (MOIN CC), Medicine, Baltimore, MD, USA. Cambridge, United Kingdom.

14Jlvmi Consulting Llc, Dousman, Wi, Usa

#See Acknowledgements for a list of all HP 13C MRI Consensus Group Members This work was supported by the ISMRM Hyperpolarized Media MR Study Group, the ISMRM Hyperpolarization Methods & Equipment Study Group, and the Hyperpolarized MRI Technology Resource Center (NIH/NIBIB grant P41EB013598).

Abstract

MRI with hyperpolarized (HP) 13C agents, also known as HP 13C MRI, can measure processes such as localized metabolism that is altered in numerous cancers, liver, heart, kidney diseases, and more. It has been translated into human studies during the past 10 years, with recent rapid growth in studies largely based on increasing availability of hyperpolarized agent preparation methods suitable for use in humans. This paper aims to capture the current successful practices for HP MRI human studies with [1-13C]pyruvate - by far the most commonly used agent, which sits at a key metabolic junction in glycolysis. The paper is divided into four major topic areas: (1) HP 13C-pyruvate preparation, (2) MRI system setup and calibrations, (3) data acquisition and image reconstruction, and (4) data analysis and quantification. In each area, we identified the key components for a successful study, summarized both published studies and current practices, and discuss evidence gaps, strengths, and limitations. This paper is the output of the “HP 13C MRI Consensus Group” as well as the ISMRM Hyperpolarized Media MR and Hyperpolarized Methods & Equipment study groups. It further aims to provide a comprehensive reference for future consensus building as the field continues to advance human studies with this metabolic imaging modality.

Keywords: Hyperpolarized MRI, metabolic imaging, carbon-13, pyruvate, dissolution dynamic

Introduction

MRI with hyperpolarized 13C agents, also known as hyperpolarized (HP) 13C MRI, has shown great potential as a novel imaging modality, particularly for its ability to probe metabolic processes in real time. The first human studies with HP [1-13C]pyruvate were performed in 2011 in prostate cancer patients (1).

Since then, there have been over 60 papers published with imaging results of human subjects from 13 different sites, with applications including prostate cancer, brain tumors, breast cancer, kidney cancer, pancreatic cancer, metastatic disease, liver disease, ischemic heart disease, diabetes and cardiomyopathies. The vast majority of these studies used [1-13C]pyruvate (1–63), where [2-13C]pyruvate (64) and 13C-urea (56) have been demonstrated too.

As clinical HP 13C MRI advances, there is a growing need to build consensus for best practices, which are critical for comparing data across sites, performing multi-site trials,deploying methods to new sites, partnering with vendors, and potentially for obtaining broader regulatory approvals.

In March 2022, we initiated an effort to build consensus within the HP 13C MRI community with this opportunity in mind, and it was greeted with strong enthusiasm. The “HP 13C MRI Consensus Group”, containing over 55 members from 27 sites, identified the area of greatest need and opportunity for consensus building to be HP [1-13C]pyruvate human

●

Pyruvate is the most mature and widely used HP agent and has the most significant translational evidence emphasizing the potential clinical impact.

●

Clinical trials, particularly multi-site trials, have the strongest need for consensus methods to ensure that data can be combined across sites. This work is a Position Paper for which the goal is to describe current successful practices and study methods for HP [1-13C]pyruvate human studies along with justification to support those practices. This is divided into four major topic areas: (1) HP 13C-pyruvate preparation, (2) MRI system setup and calibrations, (3) data acquisition and image reconstruction, and (4) data analysis and quantification (Fig. 1). The current successful practices and study methods include a literature review of published peer-reviewed journal papers showing human HP [1-13C]pyruvate study data, up to September 2022 (1–63), as well as new unpublished information from surveys of HP 13C study sites. Based on this information, we also highlight the evidence gaps, strengths, and limitations of current practices which are summarized at the end of each section.

Figure 1: Illustration of the HP 13C MRI human study process, including the 4 major areas covered in this paper: Hyperpolarized 13C-pyruvate preparation, MRI system setup and calibration, Acquisition and Reconstruction, and Data Analysis and Quantification.

Figure 2: Anatomical targets of HP [1-13C]pyruvate MRI human studies published up to September 2022.

Hyperpolarized 13C-Pyruvate Preparation

This section covers the processes for creating the HP agent, 13C pyruvate, and will include many aspects and considerations that are needed to safely and effectively prepare doses for metabolic imaging studies in human subjects. These include material, personnel, equipment and facility, fluid path preparation, quality control, and release.

It is helpful to understand that the specifications of a dose of 13C pyruvate suitable for in vivo MR HP metabolic imaging were shaped in part by early preclinical studies performed by GE HealthCare summarized in Ref. (65). In short, the safety of the two novel drug components, 13C pyruvate and the electron paramagnetic agent (EPA) AH111501, were demonstrated in those studies. The more precise formulation of the dose suitable for human use was then determined from clinical studies (66) that included two Phase 1 clinical trials in young and elderly healthy volunteers without hyperpolarization of the 13C nuclei and another Phase 1/2a dose escalation and imaging feasibility study with HP 13C pyruvate in 31 prostate cancer patients at the With the exception of the first HP 13C imaging clinical trial, which utilized a prototype device in a cleanroom (1), all HP 13C studies performed in humans to date have utilized the SPINlab polarizer (manufactured by GE HealthCare). Consequently all doses of the HP 13C pyruvate delivered by SPINlab have been produced using the “SPINlab Pharmacy Kit” that serves as the container-closure system for the various drug components (13C pyruvic acid and EPA mixture, dissolution medium, and neutralization and dilution medium) during sample polarization, dissolution and quality control (QC) processes. Thus many aspects of the HP sample preparation considerations discussed below are related to the SPINlab instrument and the consumables designed to be used with it (67).

General Considerations

While more than 860 patients or healthy subjects having been injected with HP 13C pyruvate as of January 2022 without reports of any serious adverse events (68), HP 13C pyruvate injection remains an investigational MR contrast agent and can only be administered by those with Investigational New Drug (IND) exemption from the Food and Drug Administration (FDA) in the USA, a Clinical Trial Application (CTA) in Canada, approval from National Research Ethics Committee Services in the UK, or approval from the relevant local regulatory body. Thus, methods and processes involved to produce a dose should have patient safety as the first priority. Since utilizing dissolution dynamic nuclear polarization (dissolution-DNP) for human use is still a relatively new development, there are no existing published regulatory guidelines specifically for this method.

There are two major production styles that determine how various sites approach the agent preparation. In the US, the most common approach is to rely on a sterilizing filter (“Terminal Sterilization”) to ensure sterility of the final product, akin to PET tracer production, where a starting molecule with a radioisotope is processed using various other ingredients to make the final, desired and injectable contrast agent within a necessarily short amount of time (69). For these sites, sterilization of the components and accessories upstream of this filter are not required, although many of them were manufactured and tested following Good Manufacturing Practice (GMP) or Good Laboratory Practice (GLP) requirements. The filling process is usually performed under an ISO 5 laminar flow hood, but a clean room or an isolator is not required.

This approach is typically accompanied by testing the integrity of the sterilizing filter prior to release of the dose for injection. Typically, post release endotoxin and sterility tests are performed using an aliquot reserved from each released dose.

In the UK and EU, the most common approach is to more-closely follow sterile pharmaceutical compounding guidelines (70), where all components and ingredients are required to be sterile or manufactured under GMP guidelines and are assembled and filled within a clean room environment or an isolator system (“Sterile Preparation”). Typically a batch of Pharmacy Kits for HP 13C pyruvate injection are prepared together. The sterility of the final dose is also ensured by batch validation testing, in addition to the sterility of the ingredients and the sterile compounding process. The endotoxin and sterility testing are performed for the process validation but are not performed for each injected dose.

Some institutions fill and assemble the Pharmacy Kit required for a specific study on the same day or the day prior to polarization, dissolution, and patient administration, but others have also demonstrated the feasibility of preparing a batch of kits, keeping them in a -20ºC freezer and using them over a period of a few months.

Beyond the obvious requirements that the process and the facility has to ultimately produce a dose that is safe to inject into a human, regulatory authorities will also focus on the question “Are you in control of your processes?”. To be in control of your process requires an in-depth and broad understanding of all processes involved in pre, post, and during the production process.

Personnel

It is typical and may be required to have licensed personnel involved in the production process depending on local regulations.Typically a pharmacist, radiopharmacist or other similarly qualified person (QP), in charge of the facility where the Pharmacy Kit filling and preparation is taking place, is responsible for the overall process and the release of the injectable dose.

Qualified cleanroom technicians are often involved in the Pharmacy Kit filling under the supervision of the pharmacist or QP. As is required for pharmaceutical compounding or PET tracer production, training requirements and training records for all personnel need to be maintained and available for audit by the FDA or equivalent.

Equipment And Facility

The facility and all equipment need to have standard operating procedures (SOPs) that describe how equipment is used, maintained, and calibrated to comply with relevant legislation. Currently, almost all the filling of the Pharmacy Kit takes place within a compounding laminar flow hood or isolator (typically ISO 5). At some sites, the filling is conducted within a cleanroom, while at others, it is conducted in a dedicated non-cleanroom space, reflecting differences in cleanroom approach and specifications between regulators worldwide (71). Some equipment or facilities, such as the compounding hood or cleanroom, may require external certified laboratories for testing.

Material Handling

Material handling guidelines (69,70) require SOPs detailing a system to track all of the materials involved in the HP production process for a particular patient dose, similar to current good manufacturing practice (cGMP) requirements for material handling for drug compounding. This includes acceptance standards, storage conditions, amount used in the patient dose for each ingredient and materials used in the assembly of the fluid path and Pharmacy Kit. Currently some users choose to open and inspect and sometimes modify the Pharmacy Kits upon arrival, but some users keep them in the sealed packaging until they are required for dose preparation.

Pharmacy Kit Filling And Assembling

As required by an IND or its equivalent, the preparation of the doses of HP 13C agent are detailed in the Chemistry, Manufacturing, and Control (CMC) section of an applicable regulatory submission; an example of this has been made available (72). It describes the processes of filling the Pharmacy Kit with the different components that make up the final drug product, and of assembling the final kit for either storage or immediate use in the polarizer. Special attention should be given to the laser welding process in order to satisfy installation qualification (IQ) and operational qualification (OQ). Typically, the final developed process is validated by process qualification (PQ) runs, during which 3 or more Pharmacy Kits are filled and used and the final HP 13C products are tested for endotoxin and sterility and to confirm that they meet the dose specifications for injections (usually including pyruvate concentration, residual EPA concentration, pH, liquid state polarization level and dose temperature). The data from 3 consecutive PQ runs are submitted as part of the IND submission (or its equivalent), and are often also reviewed by the Institutional Review Board (IRB) where the studies are conducted.

Quality Control And Dose Release

The quality control (QC) and dose release can be separated into two aspects: one is the QC and release of the filled Pharmacy Kit, and second is the QC and release of the HP 13C agent for injection, after polarization and dissolution. For institutions filling a batch of kits and storing them to use over a period of time, typically the batch can be released based on initial validation, environmental monitoring data from the day of kit production, and if filters are used during preparation of any of the components, filter integrity testing. But in some cases one or more kits are used for validation before the batch of kits are released for future use. For institutions that fill only the kits required for specific studies shortly before the experiment, the filled kits often do not go through separate release tests before they are used.

The quality control of the HP 13C pyruvate solution post dissolution is primarily performed to ensure that the agent meets the dose specifications (Table 1) before it is administered to the subject. These specifications target both safety (pH, residual EPA, temperature) and efficacy (pyruvate concentration, polarization, volume). Typically, the pyruvate concentration, residual EPA concentration, pH, dose temperature, dose volume, and liquid state polarization are measured by the QC accessory associated with the SPINlab polarizer. Some users perform a secondary measurement for one of the parameters, such as pH, using a different instrument or pH paper. For sites that do not go through a separate release testing process for batch filled kits, the integrity of the sterilization assurance filter, a part of the Pharmacy Kit, is typically tested as a part of the dose release. It is also common for these users to preserve an aliquot of the final HP 13C pyruvate solution for post-release endotoxin and sterility testing. This testing cannot be completed fast enough to test an individual dose prior to injection, but this is why other processes such as PQ runs and validation testing are done to minimize the chance a subject could be injected with a contaminated dose.

The Final Dose Release And Injection

should be done under the supervision of a licensed professional, based on local regulations.

Some Key Challenges

Many of the challenges associated with HP 13C pyruvate preparation can be attributed to the conditions required for the dissolution-DNP method of high magnetic field (~3-7 T) and very low temperature (~1 K) during polarization, with pressurized and superheated water necessary for the rapid dissolution event. These extreme conditions are quite challenging for the design of the container-closure and fluid path system. In particular, the cryogenic temperature in the polarizer requires special attention to any moisture or ambient (moist) air introduced into that portion of the fluid path, which can form an ice block at ~1 K. This ice can lead to flow restriction during the dissolution event and reduce the strength of the laser welded bond between the cryovial and its cap. This can ultimately produce failures in the dissolution step, including variations in final pyruvate concentration and pH that may fail to meet QC release criteria as well as fluid path ruptures that provide no available dose and result in polarizer down-time.

The polarization of the HP 13C pyruvate sample decays quickly over the span of a few minutes after dissolution, and thus the process of dissolution, QC for release, and injection should be completed as fast as possible to preserve the high polarization level achieved. Any delays in the preparation process, such as transportation time or equipment malfunction, can significantly reduce the final polarization and result in lower quality imaging data.

Current Practices

A summary of data collected from all sites performing clinical trials with HP 13C-pyruvate is shown in Fig. 3 and Table 1, including the specification of the final dose and how the quality control and release of the final dose are performed. There is a split in the Production Style, described in the General Considerations section above, with 8/13 sites using Sterile Preparation versus 5/13 using Terminal Sterilization. While many of the dose specifications show notable differences in acceptable ranges, all of these variations listed in tables have been successfully and safely been used to perform HP 13C pyruvate studies in humans. Their differences depend on the institutions’ preferences, resources and their particular regulatory situation. There is high similarity in pyruvate ranges, temperature ranges, EPA limits, and volume limits. There is modest variability in pH ranges and large variability in the endotoxin test limit. There is a 3-fold difference in acceptable polarization levels, which are measured to ensure a futile dose is not injected since the polarization is directly proportional to SNR. This reflects the decision by several sites to believe that useful data can be still be obtained with suboptimal polarizations.

Figure 3: Hyperpolarized agent preparation methods reported by sites currently performing HP

In House

Table 1: HP 13C-pyruvate preparation parameters, methods, and dose specifications used for quality control testing and release as well as validation. These were obtained from a survey of all sites performing clinical trials with HP [1-13C]pyruvate. The parameters used for product release are noted in bold text, otherwise these parameters are measured for batch validation or other QC measurements. The endotoxin and sterility testing are performed during process validation of the batch and/or post-injection, and largely depends on the agent production approach.

Summary

The overall safety record of HP 13C-pyruvate has been very strong, and the SPINlab hyperpolarizer has proven to provide high polarizations at human sized doses while meeting numerous QC and release criteria. A weakness remains the failure modes of the SPINlab Phamacy Kits (e.g. ice blocks, path ruptures), which are placed under extreme requirements particularly during dissolution. The preparation process still requires a high degree of expertise.

Therefore, there is a significant need to improve the reliability, robustness, and ease of operation for generating HP 13C-pyruvate doses for human studies. Furthermore, there is a divide between manufacturing and sterile compounding style preparation as well as other site-specific practices, resulting in variations in SOPs and justification required to relevant regulatory bodies. There have also been no comparisons between these approaches. It is also unclear what release criteria and QC parameters are truly required to ensure patient safety.

However, all of the reported methods are acceptable and approved by the appropriate regulatory authorities, and have led to the rapid expansion of successful human studies in recent years.

Mri System Setup And Calibrations

This section covers the MRI system setup, including the imaging system, RF coils, phantoms, and prescan calibration methods.

Imaging System

The main prerequisite for a given MRI scanner to be capable of supporting studies with HP 13C is its “broadband” capability to transmit and receive radiofrequency (RF) signal at the frequency of 13C, which is around 4 times lower than 1H. This does not come as a default on clinical MR devices. The transmit power of the broadband amplifier should also be sufficient to support the intended flip angle and RF pulse shape with the employed transmission RF coil(s) for 13C. Most studies to date use relatively low flip angles (< 90 degrees) for HP 13C in order to preserve polarization for time-resolved imaging. The capability to receive 13C signal on multiple channels is also desirable to increase SNR, as discussed further in the “RF coils” section.

The choice of magnetic field strength is primarily dependent on the metabolites’ frequency separation due to chemical shift dispersion and 1H imaging. High field strengths do not enhance hyperpolarized 13C signal as they do for 1H because the signal strength in a HP experiment relies on manipulating the population of quantum energy states outside of the MRI scanner.

However, the injected HP 13C-pyruvate and its metabolic products have greater frequency separation at higher fields, and it may thus be easier to separate and quantify these resonances at higher fields. This comes at the cost of a reduction in the achievable T2* and often reduced T1. As the initial polarization is independent of the imaging field strength it has been proposed that the increased T2* at 1.5T can potentially be exploited to increase SNR by adapting the acquisition bandwidth or reduce off-resonance imaging effects in cases when the decay of the transverse magnetization is dominated by T2* (73). In practice, 3T has been used in all published human 13C-pyruvate studies surveyed (Supporting Table S1), and comprises the majority of scanners currently in use for human studies (Table 3). A field strength of 3T is well-suited for 1H MRI anatomical reference and correlative imaging.

Stronger and more rapidly slewing magnetic field gradients support more rapid spatial encoding, particularly for metabolite-specific single-shot imaging using echo-planar imaging (EPI) or spiral imaging (See “Acquisition and Reconstruction”). Although the spatial resolution acquired for HP 13C imaging is typically much coarser than for 1H MRI, the factor of ~4 in gyromagnetic ratio leads to the same reduction factor in performance of the gradient system, so 13C experiments are potentially more limited by gradient hardware performance. To date, all human studies have used the commercially-available integrated gradient systems provided in clinical MRI scanners.

Optimization of scanner design has understandably focused on minimization of artifacts in 1H MRI, where devices such as room lights, the gradient amplifiers, and the motors driving the patient bed are checked to ensure that they do not produce RF interference at the 1H frequency, but artifacts may arise at other frequencies. Eddy current compensation is also not always appropriately adjusted for nuclei at other frequencies (74). In order to optimize for 13C, many sites have performed checks on phantoms for RF interference, gradient artifacts, and eddy currents (74), including the use of post-hoc gradient impulse response function characterisation and correction, and some vendors have fixed these issues as well.

Rf Coils

For HP 13C imaging studies in humans, RF coils for both 1H and 13C nuclei are needed, with 1H MRI providing an anatomical reference for registration and optional additional multiparametric MRI readouts. At the Larmor frequency of 13C nuclei, the relative contributions from coil noise compared to sample noise increase compared to 1H (73,75), although sample noise still is likely the dominant contributor for human-sized coils at 32.1MHz - the resonance frequency of 13C nuclei at 3T.

The key requirement for human 13C-pyruvate RF coils are that the coil geometry and sensitive volume must cover the volume of interest in the subject. Table 2 and Figure 4 shows coil configurations that have been used and optimized for applications in different anatomic regions.

Volume resonators are most commonly used for transmit, as they surround the subject to

Provide B1 Transmit Across The Fov (B1

+). While 1H relies on a large birdcage (“body”) coil built into the scanner, 13C transmit coils must be placed inside the bore. This takes up valuable space within the magnet, and also has led to the use of designs with relatively inhomogeneous

B1

+. Many human studies have used Helmholz pair resonators for transmit, including the “clamshell coil”, which has a notably inhomogeneous B1

+ Profile But Has Been Used Because Of

relatively easy integration into the scanner bore. B1

+ Variation Results In Variations In The Flip

angles that control the use of the hyperpolarized magnetization and creates errors in common HP metrics (9,76). The exception are head coils, where birdcage designs with highly

Homogeneous B1

+ can be placed around the head while easily fitting inside the bore. As with 1H MRI, higher SNR can typically be achieved by smaller receive coil elements, such as surface coils or phased arrays, and the majority of 13C receive coils used have layouts similar to 1H phased arrays.

RF coil quality control is important to ensure proper functioning of the coils to provide consistent imaging quality, especially with limited natural abundance 13C signal in vivo. It typically involves 1) a physical integrity check of the coil cables and connectors and 2) phantom SNR tests to check the coil’s performance and to monitor it over time (see Phantoms below). An useful reference for RF coil quality control is outlined in the MRI accreditation program of the American College of Radiology (77) and can be adapted for 13C coils.

Notably, configurations for brain and prostate studies used dual-tuned 1H/13C coil designs, which greatly simplify workflow and registration of 1H and 13C images, as no switching of coils is needed.

Table 2: RF coil configurations reported for human HP [1-13C]pyruvate studies.

Tx = Transmit

coil, RX = receive coil. The commonly used “clamshell” TX coil is a Helmholz pair design. For 1H RF configurations, all used the Body coil for TX unless otherwise noted, and “repositioned” indicates the 13C coil was removed for 1H imaging. One representative reference is listed for each configuration. The RF coil configurations reported in the reviewed papers are shown in Supporting Table S1.

Figure 4: Examples of RF coil configurations used for human HP [1-13C]pyruvate brain studies. (A,B) 13C Clamshell TX (Helmholz pair) and 2× 4-channel paddle RX arrays. (C) 13C Birdcage volume TX and 32-channel RX array (RX array slides into TX coil). (D) 13C Birdcage volume TX and 24-channel RX array, combined with a 1H 8-channel RX array. Image reproduced with permission from Ref (16).

Phantoms

Since hyperpolarized magnetization is non-renewable, phantoms containing 13C nuclei are important to: 1) test the multi-nuclear capabilities of the imaging system, including all parts of the signal excitation and receive chain; 2) perform calibration measurements before a scan with hyperpolarized nuclei; and 3) perform necessary pre-scan adjustments (see “Prescan Calibration” section). The phantoms currently in use are listed in Table 3. Their composition must provide sufficient 13C signal, with additional considerations of conductivity, stability, chemical shift(s) present, potential for dynamic imaging, and cost. The phantom geometries are typically either compact, in order to be used alongside the subject during a HP scan, or large enough to mimic the inner volume of a RF coil for system testing.

One popular compact design contains enriched 13C-urea at high concentration, typically 8 M, which provides a single resonance, placed inside a small container ~1 mL. The most common recipe mixes 13C-urea in a 90% water/10% glycerol solution, with glycerol used to increase the urea solubility and doping with a Gd-based contrast agent to shorten T1 which increases the potential SNR per unit time. For example, when Dotarem is added at a 3:1000 volume ratio the 13C-urea T1 is around 500 ms and T2 is around 100 ms. However, when testing pulse sequences influenced by T1 and T2, doping should be used carefully. This phantom is suitable for frequency calibration, transmit gain calibration, sequence testing, and as a fiducial marker when placed next to a patient. However, enriched 13C-urea has a relatively high cost compared to natural abundance compounds.

For larger volumes (>100 ml), the phantoms most often used contain undiluted ethylene glycol, glycerol, or dimethyl silicone. These compounds have sufficiently high carbon concentrations to provide sufficient 13C signal even with the 1.1% natural abundance of 13C. These larger phantoms matching the inner volume of an RF coil are useful for coil testing, including transmit

+) And Receive (B1

-) coil profile mapping, as well as to mimic acquisitions using in vivo FOV requirements. In this case, size and conductivity should match the expected subject size in order to mimic coil loading and get a realistic estimation of B1+. Large-volume natural abundance urea phantoms have also been used by some sites, but suffer from higher conductivity compared to biological tissues. Typically, it is easier to increase the conductivity and hence coil loading of the non-conductive phantom by adding NaCl to match physiological loading (16,78).

Dynamic phantoms that aim to mimic metabolite kinetics have also been developed (79–81), and have the potential to more closely mimic the HP experiment, but so far these are not widely used.

Prescan Calibration

Prior to performing an MRI acquisition, the so-called prescan procedure is used to set the shim parameters to maximize B0 homogeneity over the field of view (FOV) or a specific region of interest (ROI), the scanner center frequency (CF), the RF transmit gain, and the receiver gain.

While this calibration procedure is usually automated for 1H, the lack of sufficient natural abundance 13C signal prevents use of automated methods. (Although natural abundance 13C lipid signal has been detected, there are so far no reports on using this signal for prescan.) Table 3 shows current practices across sites.

Maximizing B0 homogeneity is independent of the nucleus and is therefore performed prior to 13C imaging using the 1H water signal and existing shimming tools, such as by a standard automated process (“Auto Shimming”) or using high order shimming routines. Similarly, the 13C CF can be calculated from the 1H CF using a predetermined scaling factor that depends on the target chemical shift (82). Another common approach used is to have a small, high-concentration 13C phantom, e.g. 8M 13C-urea, integrated in the RF coil or placed next to the scan subject (1). The reference frequency can also be based on real-time measurements after the HP injection but prior to imaging (83). Both the CF and B0 shimming are critical when using spectrally-selective RF pulses, as inmetabolite-specific imaging methods, where the desired excitation bandwidths are typically very narrow and frequency offsets can lead to a failure mode that is only apparent after injection.

The calibration of the RF transmit power is typically performed on a small, high-concentration 13C phantom placed near the region of interest during the scan or on a large 13C phantom of similar size and coil loading as the subject, prior to the subject scan. Reference power is often done by sweeping the power in a pulse-acquire sequence (53,62), or the Bloch-Siegert method (52,84). When using a small phantom, the location of the phantom, B1

+ Inhomogeneity As Well

as any shielding effects, e.g., when the phantom is integrated into a coil (1), may degrade the accuracy. Other methods include real-time Bloch-Siegert method measurements after the HP injection (83), and using the stronger natural abundance 23Na signal that is close enough to the 13C resonance frequency to be detected by 13C coils (82).

The receiver gain is predetermined, either systematically based on independent phantom measurements and assuming the dose and polarization of the HP compound is known prior to injection, or based on past HP imaging studies.

Power [Kw]

Phantom(s) - during study Phantom(s) - before study 13C Frequency

13C-bicarbonate doped with dimethyl silicone, various

Maximum Values

Table 3: Summary of the imaging systems, phantoms, and prescan procedures used at sites currently performing HP 13C-pyruvate human studies. These were obtained from a survey of all sites performing clinical trials with HP [1-13C]pyruvate. *Previously performed studies with a Siemens 3T Tim Trio. The imaging systems, phantoms, and prescan procedures reported in the reviewed papers are shown in Supporting Table S1.

Summary

Commercially available 3T MRI systems are by far the most commonly used for human HP 13C-pyruvate studies, although a systematic investigation of the impact of B0 has only recently been investigated (73). The multi-nuclear RF transmit and receive chain has proven sufficient for current acquisition strategies, although many sites have observed artifacts due to RF interference, gradient interference, and residual eddy currents when operating at the 13C frequency. A variety of 13C RF coils, tailored for numerous anatomical targets, have been successfully demonstrated, with the main limitation that most transmit coils take up a lot of additional space inside the bore and provide relatively inhomogeneous B1

+ Profiles. The

phantoms used have converged into generally 2 categories - small phantoms containing 13C-enriched compounds that can be used during the study and human-sized phantoms containing compounds with high carbon concentrations but without 13C enrichment that are used to test and calibrate the coils. There are no standardized compositions or geometry, and dynamic phantoms that recapitulate in vivo kinetics would be desirable but are still an emerging area. Prescan calibration procedures were not well defined in most publications, so we surveyed individual sites to determine current practices. Calibration procedures for the B0 field (13C CF and shimming) for most sites take advantage of 1H signal and methods, while methods

For Calibration Of B1

+ is more variable across sites, likely a reflection of remaining challenges in how to perform this calibration. Standardization of both phantoms and calibration procedures would synergistically improve the robustness and reproducibility of HP 13C studies.

Acquisition And Reconstruction

Data acquisition strategies in human HP [1-13C]pyruvate MRI studies must account for multiple chemical shifts, efficiently utilize the non-renewable HP magnetization, and acquire data quickly relative to metabolism and relaxation decay processes. These studies require spectral encoding to separate metabolites, necessitating pulse sequences that efficiently encode up to 5D data (3 spatial + 1 spectral + 1 temporal dimension). RF pulses must efficiently sample without immediately saturating the non-renewable HP magnetization, and sequences must acquire data quickly and be robust to both experimental and physiologic variation (e.g. B1

+ Inhomogeneity,

variation in perfusion) to ensure reproducibility and minimize scan-to-scan variability. This section covers current successful practices for data acquisition in human [1-13C]pyruvate studies, and accompanying 1H imaging, from different anatomic regions, including scan parameters and image reconstruction.

Acquisition And Reconstruction Methods

The acquisition methods used in human [1-13C]pyruvate studies can be classified into 3 categories: 1) MR spectroscopy or MR spectroscopic imaging (“MRS/I”), 2) chemical shift encoding methods, and 3) metabolite-specific imaging (Fig. 5).

Mrs/I Methods Specifically

resolve a spectrum that can be analyzed to extract expected as well as unexpected resonances, making this approach very robust. It was used in many initial studies (1).

Chemical Shift

encoding methods, most commonly the Iterative Decomposition of water and fat with Echo Asymmetry and Least-squares estimation (IDEAL) method, use imaging sequences acquired with multiple TEs and rely on a model-based separation of expected chemical shifts (85).

Metabolite-specific imaging methods use specialized RF pulses that are spatially and spectrally selective to excite individual metabolites which are then typically imaged with fast k-space trajectories such as echo planar imaging (EPI) or spirals (86).

Their Application To Different

organ systems is described below. The image reconstruction methods used in human [1-13C]pyruvate studies have typically been conventional methods (e.g. FFT, non-uniform FFT, or equivalent). The incorporation of accelerated imaging and advanced reconstruction methods including parallel imaging (4,57,87) and compressed sensing (7) has also been applied in human studies for improved spatial resolution, temporal resolution and coverage, but have the potential for additional artifacts as well as SNR losses due to ill-conditioning of the reconstruction (e.g. g-factor).

The Majority Of

published studies do not use accelerated imaging indicating the resolution and coverage achievable without acceleration is currently adequate for successful data collection. Performing coil combination, even with fully sampled data has also been shown to have specific challenges for HP human images: using naive sum-of-squares methods suffer from high noise amplification in the relatively low SNR regime of HP [1-13C]pyruvate (compared to 1H), motivating several HP 13C-specific methods that include data-driven coil sensitivity estimation which have shown obvious improvements over sum-of-squares (11).

More recently denoising techniques have been applied as post-processing of human HP data(41,42,44). The techniques applied are based on spatial-temporal singular value decomposition for unsupervised estimation of signal and noise components. They have shown improvements in apparent SNR in the brain and liver, while care must be taken to choose parameters such as the rank threshold to avoid oversmoothing and overfitting to the estimated signal components.

Prostate Studies

Prostate cancer was the first human application of HP [1-13C]pyruvate (1), and data was acquired with MRS/I methods: 1D dynamic MRS, single-slice 2D dynamic echo-planar spectroscopic imaging (EPSI), and single time point 3D EPSI. Advances in imaging strategies led to the development and application of new acquisition schemes, including undersampled 3D EPSI with compressed-sensing (7), model-based chemical shift encoding methods that use a priori information (47,59), and metabolite-specific EPI (10), all of which can provide volumetric whole-organ coverage and dynamic acquisitions.

The pyruvate bolus arrival in the prostate can vary by ± 10 s between patients, necessitating dynamic imaging to reliably and consistently capture the pyruvate bolus (18). For this reason, all currently ongoing studies acquire dynamic data. While MRS/I, chemical shift encoding, and metabolite-specific imaging can all achieve dynamic imaging, chemical shift encoding and metabolite-specific imaging provide greater dynamic and volumetric coverage (85). For scan prescriptions, the FOV is designed to provide full prostate coverage and typically to match the orientation of the anatomic imaging used for registration. Flip angles used in current studies are constant through time, as quantification with a variable-through-time flip scheme is highly sensitive to bolus timing (8) and errors in the RF transmit (B1 +) field (76).

Heart Studies

Data acquisition methods for 13C imaging in the heart must be designed to meet the demands of significant cardiac motion and blood flow. To cope with the periodic cardiac motion, most human heart studies to date used gating to the diastolic window, the longest cardiac cycle interval, which has reduced motion (2,22,28,30,35,36,38,45,52). The duration of the diastolic window limits the available data sampling time, making cardiac acquisitions the most time-constrained of the HP 13C MRI applications. The most common acquisition approach is metabolite-specific imaging with spiral k-space trajectories (2). Their single-shot imaging capability makes these methods particularly robust to motion effects. Furthermore, spiral k-space trajectories provide rapid k-space coverage and relatively benign flow and motion artifacts. The majority of studies have used 2D multi-slice acquisitions, but 3D encoding has also been used successfully (35).

Brain Studies

For HP 13C MRI of the human brain, the majority of studies have also used 2D (slice selective) acquisitions (10–12,14,16,28,33,40,41,44,51,53,60), with a trend toward volumetric coverage using 2D multi-slice metabolite-specific imaging. 3D metabolite-specific imaging of the whole brain, with phase encoding of the slice direction (34,57), has been shown to provide similar SNR efficiency (88) compared with multislice imaging. A number of studies have employed MRS/I (5,6,29,31–33,50,55) resulting in a spectrum from each voxel, which has the advantage of not requiring a priori information about which peaks to encode. This was important in early brain studies when it was not known which peaks would be detectable. Chemical shift encoding, using a set of images with different echo times and an iterative reconstruction of the individual resonances (i.e. the IDEAL approach (85)), has also been used (12,49,54), with the drawback that coverage in the slice direction was limited due to the time required to acquire multiple echo time images.

Abdomen And Breast Studies

The fundamental approaches to data acquisition and reconstruction in the abdomen and breast are largely similar to the aforementioned applications, but demand attention to particular challenges associated with these anatomic regions, especially relating to respiratory motion.

Although it has been shown that a basic 2D MRSI approach based on phase encoding and FID readout can be successfully applied for HP 13C imaging in breast (15) and kidney (13), major advantages in terms of spatiotemporal resolution and coverage have been realized using tailored approaches based on metabolite-specific imaging (43,62) and chemical shift encoding (43), which have facilitated multi-slice or 3D dynamic acquisitions over large FOVs in the abdomen (4,37,46).

The significant respiratory motion encountered in these regions can directly blur 13C images, and has further favored these rapid acquisition strategies. Motion also degrades B0 homogeneity, which can shift frequency-selective excitation profiles and introduce artifacts into rapid imaging readouts. This makes accurate determination of the acquisition center frequency and shimming essential in these regions which often cover large FOVs. (See “Prescan Calibration” section for more information). In some studies, breath-holding was used to minimize motion effects and enforce frame-to-frame data consistency (42). A pragmatic and reasonably effective approach for dealing with respiratory motion during 13C data acquisition is an initial breath-hold (as long as can be tolerated), followed by free-breathing (46,62).

1H Imaging

Collection of 1H imaging data is essential both for prescribing the 13C acquisition and for interpretation of the resulting 13C data. Multi-planar 1H scouts are acquired prior to 13C acquisition to enable graphical prescription of the 13C imaging region. All human HP 13C-pyruvate imaging studies acquire conventional MRI scans (e.g. T1- and T2-weighted volumes) for anatomic reference, aiming to cover at least the full 13C FOV. Acquiring these anatomic scans as close as possible to the time of 13C imaging (immediately before or after) minimizes potential misregistration between the data sets. Depending on the application, other advanced 1H sequences are also acquired (e.g. diffusion-weighted imaging for cancer imaging).

When contrast-enhanced data is acquired, it is done after 13C imaging, as paramagnetic contrast agents will accelerate 13C relaxation.

Reported Study Parameters

Figures 5 and 6, and Supporting Table S2 shows the reported acquisition study parameters for human HP [1-13C]pyruvate studies published as of September 2022. Figure 5 shows a mixture of MRS/I, metabolite-specific imaging, and chemical shift encoding methods have been successfully used, where spectroscopy-based methods have become less prevalent in recent studies. Figure 6 shows the acquisition timing, including the important start time and interval/temporal resolution, is quite variable across studies.

Figure 5: Acquisition methods used in published HP [1-13C]pyruvate human studies published up to September 2022, classified into: MR spectroscopy and spectroscopy imaging (MRS/I); chemical shift encoding methods, such as IDEAL, that use multiple TEs and model-based reconstructions; and metabolite-specific imaging methods that use spectrally-selective excitation to image a single resonance at a time.

Figure 6: Temporal acquisition characteristics reported in HP [1-13C]pyruvate human studies published up to September 2022. (a) Reported referencing of acquisition start times.

(B)

Acquisition start times reported when using dynamic imaging and when timing was reported relative to the end of the injection. (c) Temporal resolutions. “Not Applicable” indicates dynamic imaging was not used.

Summary

Three general categories of acquisition strategies have been used successfully for human HP 13C-pyruvate studies: MRS/I, model-based chemical shift encoding (e.g. IDEAL) methods, and metabolite-specific imaging methods. These have enabled successful studies in the prostate, heart, brain, abdomen, and breast. Recent studies increasingly have used the imaging-based strategies of metabolite-specific imaging and chemical shift encoding which are the fastest methods, although a heads-to–head comparison between techniques has not been performed.

Metabolite-specific imaging is quite popular because of its speed and compatibility with single-shot imaging, but is sensitive to B0 field variations and thus requires careful calibrations. Nearly all studies surveyed acquired data dynamically, allowing measurement of the bolus and metabolite kinetics. The exact timings and associated flip angles vary quite widely across reported studies, with no consensus yet as to how to choose these parameters. Image reconstruction is typically done directly using Fourier Transform methods, and accelerated imaging strategies are uncommon.

Data Analysis And Quantification

This section covers the analysis of data from human HP [1-13C]pyruvate studies, including modeling and metrics, visualization, as well as considerations for how to store data and metadata. Depending on study design, the analysis may need to give quantitative or semi-quantitative output reflecting a biological process or may just reflect a contrast between different regions of interest for quantitative evaluation.

Metrics

Figure 7: HP [1-13C]pyruvate raw data (A) have typically been quantified using four categories of metrics depending on the acquisition. Data acquired as a single time point are often quantified using normalized metabolite images or metabolite ratios (B). Dynamic data can be quantified using normalized metabolite images or metabolite ratios (B), or with metabolite timings such as time-to-peak (TTP) or pharmacokinetic (PK) models (C). The latter two require the data to be time-resolved. [1-13C]alanine and 13C-bicarbonate are analyzed similarly to [1-13C]lactate but omitted here for display.

Metabolite images are commonly used as summary metrics for HP MRI data, often including some form of normalization as well as summed over time as an area under the time curve (AUC) (17). These are analogous to the visual evaluation that is most used for routine clinical work (89,90). In these metabolite images, we expect that the [1-13C]pyruvate AUC signal is predominantly weighted towards perfusion and uptake, while [1-13C]lactate, [1-13C]alanine and 13C-bicarbonate AUCs represent metabolic conversion. The strength of this approach lies in its simplicity and relatively few underlying assumptions. Limitations to the use of single-metabolite images or AUCs include sensitivity to inhomogeneous coil profiles (57,87,91), the acquisition strategy and acquisition parameters, pyruvate polarization and concentration level, and signal relaxation rates (92). Further, the reader must be careful to interpret all the images in conjunction to better understand the underlying biology; for example, increased [1-13C]lactate in the presence of decreased [1-13C]pyruvate delivery can have a very different meaning compared to increased [1-13C]lactate with increased [1-13C]pyruvate delivery.

In an attempt to address variations in coil sensitivity, polarization level, and pyruvate delivery, AUC images are often computed by normalizing to a specified parameter, such as the maximum pyruvate or average lactate signals, or presented as a ratio such as lactate/pyruvate or divided by “total Carbon” - the sum total of HP 13C signal observed across all metabolites. The AUC ratios between metabolites and pyruvate are proportional to the corresponding forward kinetic rates (81,93), but are not directly comparable to rate constants when magnetization loss rates (e.g. relaxation and losses due to signal excitation) differ between studies. Similarly, the ratios between the produced metabolites (e.g. bicarbonate/lactate) can reflect the balance between downstream metabolic pathways (12,55). Care must be taken to consider how AUC images are calculated and normalized before comparing values between studies.

To further quantify the interpretation, pharmacokinetic (PK) modeling approaches were developed to compute the apparent kinetics of pyruvate-to-metabolite exchange (92,94–99). These yield semi-quantitative to quantitative apparent rate constants, given in s-1. Some models require a vascular input function, while others avoid this requirement (95). PK models can explicitly account for acquisition-specific details such as excitation angle and repetition time, and thus may reduce the effects of these details on quantification. An input-less model, provided in the Hyperpolarized-MRI-Toolbox (https://github.com/LarsonLab/hyperpolarized-mri-toolbox) (100) and thus frequently employed for human data, has been shown to fit well and robustly to prostate and brain data (8,20). PK models are quantitative in nature, arguably provide more relevant biological information (8,20), and appear to be reproducible across sites (51). However, rate constants derived from PK models are still apparent rates, and likely do not reflect a single biological characteristic.

Some additional considerations include whether complex or magnitude data is used, as the noise behaviors will impact the analysis differently. Additionally, cut-off thresholds or other criteria may be used to identify and avoid voxels with insufficient SNR before analysis to improve robustness (20,41).

Regardless of the analysis approach, the underlying biology is not always clearly represented by the data; instead, the metrics may be influenced by perfusion, barrier permeability, intercellular shuttles, enzyme activities, co-substrate concentrations, or combinations thereof, depending on the organ and disease of interest (19,43,94,101–103). This may be addressed by incorporating complementary information. As an example, HP 13C pyruvate data is influenced by perfusion, and thus addition of perfusion MRI could be important for interpretation (98,104,105).

All the methods outlined above have been explored in clinical studies, described in Supporting Table 3 and summarized in Figure 8. As of September 2022, approximately 52% of studies involving human subjects report rate constants derived from a PK model with a few different models reported. A nearly equal fraction (51%) of the studies report AUC ratio values.

Approximately 66% of these studies report metabolite-specific images or AUC values. About 40% report SNR values; this metric is particularly frequent in manuscripts that describe technical developments for clinical HP MRI. Approximately 16% of these studies summarize model-free metrics, and 10% report measurements from a single timepoint. Most studies report a combination of quantities.

Figure 8: Reported metrics used for analysis in HP [1-13C]pyruvate human studies published up to September 2022.

Visualization

A wide variety of approaches have been used for visualizing data from human HP 13C-MRI studies. The challenges and practical considerations are: 1) choosing the appropriate metrics to display, 2) how to encode the parameters (e.g. the colormap), and 3) choosing how to provide anatomical context and other multi-parametric data. The choice of visualization also depends on the goal which could be for diagnostic interpretation, but also quality control, reproducibility among readers and publication.

Metrics

The choice of HP 13C metrics is described in detail above. At this stage in HP 13C development where there is no standardized metric, often a combination of metabolite images and ratios or PK model parameters are shown.

Parameter Encoding

The mapping function chosen should provide an adequate, often quantitative, impression of the parameter mapped. There is a consensus in the visualization field that perceptually uniform maps are best suited to visualize continuous parameters, like the greyscale typically used by radiologists as well as other monochrome (black to blue) and color ranges (fire-type, rainbow-type) (106,107). Multi-color heatmaps have been the most frequently employed method for HP 13C data, while greyscale has infrequently been used but it ensures there is no coloring-based bias as well as facilitating later reuse (Fig. 9a). Among the color schemes employed in the clinical HP 13C literature, fire-type scheme seems to be the most common [similar to “Plasma” or “Inferno” in matplotlib.org]. Next most commonly employed is the rainbow-type scheme [similar to “Rainbow” in matplotlib.org].

Anatomical Context

HP MRI faces the challenge that it does not necessarily depict the anatomical features, similar to PET, and thus requires an anatomical reference. Most often, a grayscale anatomical image is overlaid with a HP colormap (Fig. 9c,d). This approach is very intuitive, but can skew perception as the grey-scale anatomical reference may affect the brightness of the HP data (e.g. signal in the skull). This bias does not occur when showing adjacent maps (Fig. 9a, b). Here, anatomical outlines may help to provide reference (Fig. 9b).

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.