TRAINING-FREE ROBUST MULTIMODAL LEARNING VIA SAMPLE-WISE JACOBIAN REGULARIZATION
Year: 2021

Abstract:
Multimodal fusion emerges as an appealing technique to improve model performances on many tasks. Nevertheless, the robustness of such fusion methods is rarely involved in the present literature. In this paper, we are the first to propose a training-free robust late-fusion method by exploiting conditional independence assumption and Jacobian regularization. Our key is to minimize the Frobenius norm of a Jacobian matrix, where the resulting optimization problem is relaxed to a tractable Sylvester equation. Furthermore, we provide a theoretical error bound of our method and some insights about the function of the extra modality. Several numerical experiments on AV-MNIST, RAVDESS, and VGGsound demonstrate the efficacy of our method under both adversarial attacks and random corruptions.

heading: 1 INTRODUCTION
text: Deep fusion models have recently drawn great attention of researchers in the context of multimodal learning (Vielzeuf et al., 2018; Baltrušaitis et al., 2018; Pérez-Rúa et al., 2019; Wang et al., 2020; Xue et al., 2021) as it provides an easy way to increase model accuracy and robustness. For instance, RGB cameras and LiDARs are usually deployed simultaneously on an autonomous vehicle, and the resulting RGB images and point clouds are referred to as two modalities, respectively. When RGB images are blurry at night, point clouds could provide complementary information and help to make decisions in vision tasks (Kim & Ghosh, 2019). Over the past few years, numerous multimodal fusion methods have been proposed at different levels: early-, middle-, and late-fusion (Chen et al., 2021). In early-fusion, input feature vectors from different modalities are concatenated and fed into one single deep neural network (DNN), while in middle-fusion, they go into DNNs independently and exchange information in feature space. Unlike the previous two cases, late-fusion is realized by merging distinct DNNs at their output layers via concatenation, element-wise summation, etc.
These three levels of fusion possess different pros and cons. For instance, late-fusion, the primary concern of our paper, is (i) privacy-friendly and (ii) convenient to deploy. Specifically, assume that a hospital wants to have an AI agent to judge whether a patient has a certain disease or not (Sun et al., 2020). It has to divide the complete training feature (e.g., medical records, X-ray images) of every patient and deliver them to different AI companies, otherwise, the patients’ identities will be exposed and their privacy are unprotected. This, in turn, directly rules out the possibility of applying early- or middle-fusion methods. On the other hand, the hospital could still exploit latefusion technique to generate the ultimate AI agent after several unimodal DNNs are trained by AI companies. Moreover, unlike early- or middle-fusion, many late-fusion methods could tolerate missing modality information (i.e., no need for paired data) and thus are convenient to deploy.
Although late-fusion is a mature topic in the literature, its performance under adversarial attacks (Madry et al., 2018; Tsipras et al., 2019) and random corruptions (Zheng et al., 2016; Kim & Ghosh, 2019) is rather under-explored. In this paper, we address the problem of robust late-fusion by utilizing Jacobian regularization (Varga et al., 2017; Jakubovitz & Giryes, 2018; Hoffman et al., 2019; Chan et al., 2019) and conditional independence assumption (Sun et al., 2020). The key is to minimize the Frobenius norm of a Jacobian matrix so that the multimodal prediction is stabilized (see Figure 1). Our main contributions are as follows:
• To the best of our knowledge, we are the first to propose a training-free robust late-fusion method. The involving optimization problem is relaxed to a Sylvester equation (Jameson, 1968), and the solution is obtained with only a little computational overhead.
• We provide a theoretical error bound of our proposed robust late-fusion method and an illustrative explanation about the function of the extra modality via the TwoMoon example.
• Thorough numerical experiments demonstrate that our method outperforms other latefusion methods and is capable to handle both adversarial attacks and random corruptions.

heading: 2 PRELIMINARY
text: Network Robustness To verify the network robustness, two major kinds of perturbations are used, which in our paper we referred to as (i) adversarial attacks such as FGSM, PGD, or CW attack (Goodfellow et al., 2014; Madry et al., 2018; Carlini & Wagner, 2017) and (ii) random corruptions such as Gaussian noise, missing entries or illumination change (Zheng et al., 2016; Kim & Ghosh, 2019). Correspondingly, many methods have been proposed to offset the negative effects of such perturbations. Adversarial training based on projected gradient descent (Madry et al., 2018) is one strong mechanism to defend against adversarial attacks, and the recent Free-m method (Shafahi et al., 2019) is proposed as its fast variant. Besides adversarial training, several regularization techniques are also proven to have such capability, such as Mixup (Zhang et al., 2020), Jacobian regularization (Jakubovitz & Giryes, 2018). Alternatively, with regard to random corruptions, Mixup and Jacobian regularization are also effective in this case (Zhang et al., 2018; Hoffman et al., 2019). Another powerful approach is stability training (Zheng et al., 2016), where it introduces an additional KL divergence term into the conventional classification loss so that the trained DNNs are stabilized.
Multimodal Learning DNNs trained by fusing data from different modalities have outperformed their unimodal counterparts in various applications, such as object detection (Chen et al., 2021; Kim & Ghosh, 2019), semantic segmentation (Chen et al., 2020b; Feng et al., 2020), audio recognition (Gemmeke et al., 2017; Chen et al., 2020a). Based on where the information is exchanged between different modalities, multimodal fusion methods could be classified into three kinds: (i) early-fusion (Wagner et al., 2016; Chen et al., 2021), (ii) middle-fusion (Kim & Ghosh, 2019; Wang et al., 2020), and (iii) late-fusion (Chen et al., 2021; Liu et al., 2021). For instance, if the information is fused at the end of DNNs (e.g., Figure 1), such a method belongs to the late-fusion category.
Although vast efforts have been put into exploiting multimodal fusion methods to improve DNNs’ accuracy on specific learning tasks, few works have explored network robustness in the multimodal context. Specifically, Mees et al. (2016), Valada et al. (2017), and Kim et al. (2018) exploited gating networks to deal with random corruptions, adverse or changing environments. Afterwards, Kim & Ghosh (2019) proposed a surrogate minimization scheme and a latent ensemble layer to handle single-source corruptions. However, all the aforementioned methods belong to middle-fusion and only random corruptions are considered. On the other hand, we focus on another important scenario: late-fusion, and besides corruptions, we further take adversarial attacks into account. To the best of our knowledge, robust late-fusion is un-explored in the previous literature.

heading: 3 JACOBIAN REGULARIZATION IN TEST TIME
text: Consider a supervised K-class classification problem in the multimodal context. Suppose that features from two distinct modalities A (e.g., audio) and B (e.g., video) are provided in the form of DA = {(xiA,yi)}Ni=1 and DB = {(xiB ,yi)}Ni=1, where {xA,xB} represents input features and y represents the true label. We train two unimodal networks separately, each corresponding to one modality. Given a specific input feature xA, the first unimodal network calculate the class prediction pA ∈ RK by:
pA = σ1(zA) = σ1(WAhA + bA) = σ1(WAfA(xA) + bA) (1)
where σ1(·) represents the Softmax function, zA ∈ RK represents the raw logit, and hA = fA(xA) ∈ RH represents the feature being fed into the last layer 1. Here WA ∈ RK×H and bA ∈ RK are the learnable weight and bias of the last linear layer, respectively. Similarly, the second unimodal network provides the class prediction pB = σ1(zB) = σ1(WBhB + bB) = σ1(WBfB(xB) + bB).
Based on the conditional independence assumption (Kong & Schoenebeck, 2018), the basic statistical late-fusion method generates the final class prediction as (Chen et al., 2021):
p = σ2( pA pB
freq ) = σ2( σ1(zA) σ1(zB) freq ) (2)
where the division is performed in an element-wise manner, represents the element-wise product, and σ2(·) represents a linear normalization enforcing the summation of elements equal to one. Here freq ∈ RK contains the occurring frequencies of each class, calculated from the training dataset. Our proposed approach builds upon (2). Specifically, we consider adding twoK×K weight matrices {Wa, Wb} ahead of {zA, zB}, respectively, right before they get activated. Consequently, the final multimodal prediction is re-calibrated as:
p′ = σ2( σ1(z
′ A) σ1(z′B)
freq ) = σ2( σ1(WazA) σ1(WbzB) freq ) (3)
Suppose that the data provided from modality A are perturbed at the input or feature level, while the data from modality B are clean. For matrix Wb, we could simply set it to an identity matrix, implying that we didn’t invoke the robust add-on for the second modality. To determine the value of Wa, we first calculate the derivative of p′ with respect to hA (see Appendix A.1):
∂p′ ∂hA = J′WaWA = [p ′p′,T − Diag(p′)]WaWA (4)
Then we minimize the following regularized Jacobian loss with respect to Wa:
min Wa L = min Wa
(1− γ)||J′WaWA||2F + γ||Wa − I||2F (5)
where 0 < γ < 1 is a tunable hyper-parameter. Minimizing the first term in the loss could make the change of p′ limited to some extent when hA is perturbed, while the second term in the loss guarantees numerical stability and the prediction in the perturbed case won’t get too far from that of the clean case. For a specific multimodal input {xA,xB}, once Wa is determined, so are p′ and J′ via (3) and (4), respectively. Thus, (5) is well-determined and non-linear with respect to Wa.
We propose a heuristic iterative method making the above optimization problem tractable in Algorithm 1. Our key is to decouple J′ from Wa. Namely, in step 5 of Algorithm 1, all terms are known except Wa, and thus the relaxed loss is convex. After writing out all terms of ∂L(t)/∂Wa, we observe that it is a Sylvester equation (Jameson, 1968). It has a unique solution and the run time is as large as inverting a K ×K matrix and hence affordable. See Appendix A.2 for details. Remarks First, in our implementation we find that one iteration could already yield a sufficiently accurate result, thus all our numerical results are reported with tmax = 1 (i.e., p′ = p(1)). Second, if we know data from modality B are also perturbed, we could solve Wb in a similar manner. Third,
1(1) holds because the last layer is usually implemented as a fully connected (FC) layer with Softmax activation in classification tasks.
Algorithm 1 Iteratively solving regularized Jacobian loss 1: Given one specific input {xA,xB}. Initialize iteration index t = 0. 2: Perform one forward pass, yielding the initial class prediction p(0). 3: while t < tmax do 4: Calculate J(t) = p(t)p(t),T − Diag(p(t)). 5: Minimize L(t) = (1− γ)||J(t)WaWA||2F + γ||Wa − I||2F with respect to Wa. 6: Calculate p(t+1) based on Eq (3) with the optimal Wa. 7: Update t = t+ 1. 8: end while 9: Return p′ = p(t).
notice that our approach is invoked merely during inference (i.e., test time), but not the train time. Furthermore, we demonstrate our approach in the context of two modalities, while obviously, it could be equally applied to many modalities. Finally, with a moderate assumption, we can further prove that when the input is perturbed, the change of final prediction enhanced by our method is limited to a certain amount. This has been summarized in Theorem 1 (see Appendix A.3). Some immediate corollary include: (i) when ∼ N(0,Σ), the bound is simplified toE[||p′,noise−p′||] ≤ l ( γK2(1−γ) ) 1/2 Tr[Σ], (ii) when the L2 norm of is constrained smaller than δ (usually assumed in adversarial attacks), the bound is simplified to ||p′,noise − p′|| ≤ l δ ( γK2(1−γ) ) 1/2.
Theorem 1 If fA is l-Lipschitz continuous and xA is perturbed by : xnoiseA = xA + , then the Euclidean norm of our final prediction (i.e., ||p′,noise − p′||) at most changes l √ γK
2(1−γ) || ||.

heading: 4 NECESSITY OF THE EXTRA MODALITY: BIASING EFFECT
text: In this sub-section, we explain the principles behind our approach. Let us move one step back and consider the unimodal case. In what follows, we will demonstrate that our approach might not work well in the unimodal context, which in turn justifies the necessity of the extra modality.
At first glance, our approach seems to be equally applicable in the unimodal context. Namely, suppose that only one unimodal network p = σ1(z) = σ1(Wh + b) is available, in test time and for each specific input x, we add a weight matrix Wx before the raw logit z being fed into the Softmax activation, re-calibrating the final unimodal prediction as: p′ = σ1(z′) = σ1(Wxz). Here Wx can be solved by analogy to (5). However, we observe that the introduction of Wx usually won’t change the final prediction.
TwoMoon example For an intuitive understanding, we consider the TwoMoon example used by Xue et al. (2021). In this example, all data points are scattered in a 2D space. The data located at the upper and lower leaf have true labels 0 and 1, and are colored by red and blue, respectively.
In the unimodal case, we take both horizontal and vertical coordinates as input and train a neural network with three FC layers. As shown in Figure 2 (a), the network perfectly fits to the clean data and achieves 97.75% accuracy on the clean test data. In the remaining three figures, we evaluate the trained network on noisy test data. Specifically, we deliberately choose γ = 1.0 in Figure 2 (b), so that our approach is actually not invoked (since the solved Wx equals I). In this case, the accuracy drops to 97.00% on noisy data, while the heatmap in the background doesn’t change at all compared to Figure 2 (a). In Figure 2 (c) and (d), we choose γ = 0.5 and γ = 0.1, respectively. We observe that even though the color of the heatmap surely becomes lighter as expected, the test accuracies of both cases are still 97.00%. More importantly, we find that this phenomenon is not coincident and that the prediction doesn’t change for any input after adding Wx in unimodal binary classification. See Appendix A.4 for rigorous mathematical proof. Moreover, in a K-class (K > 2) classification problem, the final prediction might change if γ is sufficiently small. For instance, given z = [1, 0, 2]T and W = I, if we choose γ = 0.5, then the final prediction will change from p = [0.245, 0.090, 0.665]T to p′ = [0.270, 0.096, 0.635]T , where the last entry remains to be the largest. However, if we choose γ = 0.01, then the final prediction will become p′ = [0.391, 0.219, 0.390]T , and now the first entry is the largest. See Appendix A.5 for a theoretical bound of γ in this highdimensional case.
Now we turn to the multimodal case. We treat the horizontal coordinates as one modality, and the vertical coordinates as the second modality. Two neural networks are trained, each corresponding to one modality. Then we fuse them based on the aforementioned statistical fusion method. As shown in Figure 3, in the multimodal context, when our method is invoked (γ = 0.5 or 0.1), the color of the heatmap becomes lighter and the test accuracies on noisy data all increase compared to the trivial statistical fusion (i.e., when γ = 1.0). On the other hand, the test accuracies on clean data almost remain still (or slightly increase alongside γ ↓).
Lemma 1 Statistical fusion is equivalent to concatenating the raw logits:
σ2( pA pB
freq ) = σ1(zA + zB − ln freq) = σ1([ I I ] [ zA zB ] − ln freq)
Biasing effect of the extra modality As we have seen in the TwoMoon example, our proposed method behaves differently in the unimodal and multimodal context. To get a better understanding of this subtle difference, we first summarize an equivalent form of statistical fusion in Lemma 1 (see Appendix A.1). Now assume that the sample-wise Wa is added, then the final multimodal prediction becomes σ1(WazA + zB − ln freq). Alternatively, in the unimodal context, our method adds a sample-wise Wx and the final prediction becomes σ1(Wxz). Comparing these two expressions,
we observe that the introduced Wa and Wx occur in the perturbed modality. However, in the multimodal case, the extra modality acts as a bias term inside the Softmax function σ1 and plays a key role.

heading: 5 EXPERIMENTAL RESULTS

heading: 5.1 AV-MNIST
text: AV-MNIST (Pérez-Rúa et al., 2019; Vielzeuf et al., 2018) is a novel dataset created by pairing audio features to the original MNIST images. The first modality corresponds to MNIST images with size of 28 × 28 and 75% energy removed by principal component analysis (Pérez-Rúa et al., 2019). The second modality is made up of spectrograms with size of 112 × 112. These spectrograms are extracted from audio samples obtain by merging pronounced digits and random natural noise. Following Vielzeuf et al. (2018), we use LeNet5 (LeCun et al., 1989) and a 6-layer CNN to process image and audio modalities, respectively. To thoroughly verify the efficiency of our method, we design various types of perturbations: (i) random corruptions including Gaussian noise (ω0), missing entries (ω1), and biased noise (ω2, ω3), and (ii) adversarial attacks including FGSM (ω4) and PGD attack (ω5, ω6). See Appendix B.1 for detailed definitions.
Baselines We limit our comparison to the range of late-fusion methods. To verify our method, we consider several methods which could improve network robustness, including (i) regular training, (ii) adversarial training (advT) by Madry et al. (2018), (iii) free-m training (freT) by Shafahi et al. (2019), (iv) stability training (stabT) by Zheng et al. (2016), and (v) Mixup (mixT) by Zhang et al. (2018). To the best of our knowledge, no previous robust late-fusion methods exist. Thus for comparison purpose, we slightly modify the confidence-based weighted summation in Yu et al. (2021) and adopt it here, referred to as mean fusion with confidence-based weighted summation. Combinations of different fusion methods and robust add-ons provide us with a few baselines. We emphasize that in our experiments, the train data are always free of noise. Following Kim & Ghosh (2019), our experiments are mainly conducted on the case when one modality is perturbed.
Results Table 1 reports the test accuracies of several models on AV-MNIST. Six unimodal networks are trained on clean data. Then we fuse them with different fusion methods and different robust add-ons. We only invoke our Jacobian regularization for the perturbed audio modality, while keeping the image modality unchanged. First, in the case of mean fusion, confidence-based weighted summation doesn’t improve robustness. In the case of statistical fusion, we observe that our proposed Jacobian regularization method not only boosts the accuracy on noisy data but also on clean data. For instance, when our method is invoked on the combination of model-0 and model-1, the test accuracy on clean data increases from 93.6% to 94.7%, and the test accuracy on data with Gaussian noise (e.g., ω0 = 1.0) increases from 80.4% to 84.5%. Similar phenomena can be observed when our method is invoked on other combinations. Moreover, the number of green arrows is much larger than that of red arrows implying that our method is compatible with different robust techniques and applicable under different noises/attacks.
Comparing the last row in the second sub-section and the third row in the third sub-section, we find that with our proposed method enabled, even if the unimodal backbone model-1 is worse than model-2, the final multimodal accuracies (0, 1) with our method invoked are even larger than (0, 2) without our method invoked. Similar phenomena could usually be observed in other sub-sections. We argue that in this example, statistically fusing two regular unimodal baselines with our Jacobian regularization invoked surpasses robust unimodal baselines and their trivial fusion. Moreover, the highest accuracy in each column has been bolded. It is clear that the multimodal network with our method invoked usually achieves the best performance over all others.
Figure 4 further plots model accuracy versus magnitude of noise/attack. It displays a consistent trend that our method works well regardless of the magnitude of noise/attack. Furthermore, the larger the noise/attack is, the better our method performs. More numerical results (such as the impact of hyper-parameters, applying perturbation on the image modality) are shown in Appendix B.2.

heading: 5.2 EMOTION RECOGNITION ON RAVDESS
text: We consider emotion recognition on RAVDESS (Livingstone & Russo, 2018; Xue et al., 2021). We use a similar network structure for both the image and audio modality: Three convolution layers, two max pooling layers, and three FC layers connect in sequential. Similar baselines and experimental settings are adopted as those in AV-MNIST. Results are reported in Table 2.
Results Still, in the case of mean fusion, confidence-based weighted mean fusion usually doesn’t work except in few cases, while our proposed Jacobian regularization exhibits good compatibility with different trained models under various perturbations. The bold results imply that the multimodal network with our Jacobian regularization invoked outperforms all other models except in FGSM attack (ω4 = 0.03). Due to page limits, extra results are presented in Appendix B.3.

heading: 5.3 VGGSOUND
text: We consider classification task on a real-world data set VGGSound (Chen et al., 2020a), where two modalities audio and video are available. To construct an affordable problem, our classification task only focuses on a subset of all classes. Namely, we randomly choose 100 classes and our classification problem is constrained on them. Consequently, there are 53622 audio-video pairs for training and 4706 audio-video pairs for testing. For the audio modality, we apply a short-time Fourier transform on each raw waveform, generating a 313 × 513 spectrogram. We take a 2D ResNet-18 (He et al., 2016) to process these spectrograms. For the video modality, we evenly sample 32 frames from each 10-second video resulting input feature size of 32 × 256 × 256. We take a ResNet-18 network to process the video input where 3D convolutions have been adopted to replace 2D convolutions.
Results We present the results of mean fusion and statistical fusion with or without robust add-on in Table 3. Unlike the previous experiments, here we make audio modality clean and assume corruptions/attacks on the video modality. The results demonstrate that our method outperforms counterpart statistical fusion without Jacobian regularization in most cases. Furthermore, the bold results imply that our Jacobian regularization could enhance model robustness under various corruptions or attacks. More results (e.g., unknown corrupted modlity) are presented in Appendix B.4

heading: 6 DISCUSSION AND FUTURE WORK
text: In this paper, we propose a training-free robust late-fusion method. Intuitively, our method designs a filter for each sample during inference (i.e., test time). Such a filter is implemented by Jacobian regularization, and after a sample goes through it, the change of final multimodal prediction is limited to some amount under an input perturbation. The error bound analysis and series of experiments justify the efficacy of our method both theoretically and numerically. The Twomoon example explains the biasing effect of the extra modality, rooted for the difference between unimodal and multimodal.
Our method opens up other directions for further exploration. In the unimodal context, we understand that directly applying our method usually doesn’t change the prediction. Thus, would it be capable to adjust the confidence of a DNN? Another possibility is that besides Wx, we deliberately add another bias term bx and optimize both for unimodal robustness. Alternatively, in the multimodal context, we might consider minimizing the derivative of the final prediction directly to the input feature. Moreover, it is educational to compare it with consistency regularization. Last but not least, statistical fusion implicitly uses the assumption that train and test data come from the same distribution, so that we could calculate freq based on the train dataset and reuse it in inference time. When domain shift presents, this doesn’t work and a straightforward remedy might be to make freq a learnable parameter.

heading: A MISSING DEDUCTION

heading: A.1 EQUIVALENT FORM AND DERIVATIVE OF STATISTICAL FUSION
text: In this subsection, we first prove Lemma 1. For conciseness, we denote the k-th element of pA as pA,k. Similar notations are adopted for pB , zA, zB , and ln freq. The k-th element of the left-hand side (LHS) in Lemma 1 can be simplified as:
pA,kpB,k/freqk∑K j=1 pA,jpB,j/freqj =
exp(zA,k) TA exp(zB,k) TB /freqk∑K j=1 exp(zA,j) TA exp(zB,k) TB /freqj
= exp(zA,k) exp(zB,k)/freqk∑K j=1 exp(zA,j) exp(zB,k)/freqj = exp(zA,k + zB,k − ln freqk)∑K j=1 exp(zA,j + zB,j − ln freqj)
(6)
which is exactly the k-th element of the right-hand side (RHS) in Lemma 1 by definition. Note that in the first line of (6), we use the definition pA = σ1(zA), pB = σ1(zB), and have introduced two normalization constants TA = || exp(zA)||1 and TB = || exp(zB)||1. By analogy to the above deduction, an equivalent form of (3) is p′ = σ1(WazA+WbzB−ln freq). Based on the chain rule and zA = WAhA + bA, we have:
∂p′ ∂hA = ∂p′ ∂zA ∂zA ∂hA = J′WaWA = [p ′p′,T − Diag(p′)]WaWA (7)
where we have used the Jacobian matrix of Softmax function to calculate J′ (Bishop, 2006). We emphasize that J′ is symmetric.

heading: A.2 SOLUTION OF REGULARIZED JACOBIAN LOSS
text: Because Frobenius norm is convex and a weighted summation with positive coefficients preserves convexity, the relaxed (5) shown in Step 5 of Algorithm 1 is a convex function with respect to Wa. This implies that any local minimum is also the global minimum. We calculate the derivative of the relaxed (5) with respect to Wa and set it to zero, yielding:
( 1
γ − 1)J(t),2Wa(WAWTA) + Wa = I
⇒ AWa + WaB = B (8)
where for simplicity we have defined A = (1/γ − 1)J(t),2 ∈ RK×K and B = (WAWTA)−1 ∈ RK×K . (8) is known as Sylvester equation2, and it has been proven that if A and −B do not share any eigenvalues, then the Sylvester equation has an unique solution (Jameson, 1968; Horn & Johnson, 2012; Bartels & Stewart, 1972). In our case, since A is positive semi-definite (PSD) and −B is negative semi-definite, the only possible overlapping eigenvalue is 0. On one hand, 0 is surely an eigenvalue of A since Ju = 0u, where u is a column vector with all elements equal to 1. On the other hand, it is argued that a well-generalized neural network normally doesn’t have singular weight matrices (Martin & Mahoney, 2021). Thus, 0 is almost surely not the overlapping eigenvalue and a unique solution of (8) is guaranteed. As demonstrated by Jameson (1968), the matrix equation (8) can be converted to invert a K × K matrix. As in a classification task, K won’t be too large, the run time of solving Sylvester equation is affordable in this context. Furthermore, as shown in the above equation, when γ → 0, the matrix A will tend to infinity, leading to an ill-conditioned matrix equation. This, in turn, implies that the second term in the loss shown in Eq (5) guarantees numerical stability.
2A general Sylvester equation has the form: AX + XB = C, where {A,B,C} are all known and X is required to be solved.

heading: A.3 PROOF OF THEOREM 1
text: In the last iteration t = tmax →∞, noticing Wa satisfies (8), we have:
||J(tmax)WaWA||2F = Tr[J(tmax)WaWAWTAWTa J(tmax)] = Tr[J(tmax),2Wa(WAWTA)W T a ]
= ( 1
γ − 1)−1Tr[AWaB−1WTa ]
= ( 1
γ − 1)−1Tr[(B−WaB)B−1WTa ]
= ( 1
γ − 1)−1Tr[Wa −WaWTa ]
(9)
and: ||Wa − I||2F = Tr[WaWTa − 2Wa + I] (10)
Thus we have:
||J(tmax)WaWA||2F ≤ 1
1− γ [(1− γ)||J(tmax)WaWA||2F +
γ 2 ||Wa − I||2F ]
= 1
1− γ
{ γTr[Wa −WaWTa ] + γ
2 Tr[WaWTa − 2WTa + I] } = γ
2(1− γ) Tr[I−WaWTa ]
≤ γ 2(1− γ) Tr[I] = γK 2(1− γ)
(11)
where we have taken advantage of the fact that WaWTa is a positive semi-definite matrix, and thus its trace is no less than zero. Based on the first-order Taylor expansion, we have:
p′,noise − p′ ≈ ∂p ′
∂hA (hnoiseA − hA)
= ∂p(tmax+1)
∂hA (hnoiseA − hA)
≈ J(tmax)WaWA(hnoiseA − hA)
(12)
Thus we have: ||p′,noise − p′|| ≤ ||J(tmax)WaWA||F · ||hnoiseA − hA||
≤
√ γK
2(1− γ) · l||xnoiseA − xA||
= l
√ γK
2(1− γ) || ||
(13)
where we have utilized the l-Lipschitz continuity of fA.

heading: A.4 UNIMODAL 2D CASE
text: Consider an illustrating example when K = 2 (a.k.a. unimodal binary classification). In this case, the final prediction is two dimensional: p = [p1 p2]T . We denote the raw logit after adding the Wx matrix as z′ = Wxz, and its entries are denoted as z′1 and z ′ 2. Let us begin by explicitly writing out the first line of (8). After some simplification, we obtain:[ κ −κ −κ κ ] [ ρ1 ρ2 ρ3 ρ4 ] [ a b b c ] + [ ρ1 ρ2 ρ3 ρ4 ] = I (14)
where κ = 1/γ − 1. Here we have explicitly written out the elements of WWT and merged some terms of J′,2 into a = p1p2(w211 + w 2 12), b = p1p2(w21w11 + w12w22), and c = p1p2(w 2 21 + w 2 22).
The unknowns occur in Wx:
Wx = [ ρ1 ρ2 ρ3 ρ4 ] (15)
After some algebra, we can solve all unknowns via the following block matrix inversion: ρ1ρ2ρ3 ρ4  = [ G HH G ]−1 ·  100 1  (16) where for simplicity, we have defined:
G = [ 1 + aκ bκ bκ 1 + cκ ] , H = [ −aκ −bκ −bκ −cκ ] (17)
Noticing that z′ = Wxz, we can prove that: (z′1 − z′2) · (z1 − z2) = zT [ ρ1 − ρ3 −(ρ1 − ρ3) ρ2 − ρ4 −(ρ2 − ρ4) ] z (18)
It is easy to calculate the eigenvalue of the matrix in the middle is 0 and ρ1 + ρ4 − ρ2 − ρ3. In other words, if we could prove ρ1 + ρ4 − ρ2 − ρ3 is no less than zero, then the matrix in the middle is positive semi-definite, which implies order preserving property holds: If z1 > z2, then z′1 > z ′ 2 and vice versa. This order preserving property is equivalent to say that adding the Wx matrix makes no sense in unimodal binary classification problem, since the predication won’t change for any input.
Now let us prove ρ1 + ρ4 − ρ2 − ρ3 ≥ 0. With a lot of algebra, we obtain:
ρ1 + ρ4 − ρ2 − ρ3 = [ 1 −1 −1 1 ]  ρ1ρ2ρ3 ρ4 
= [ 1 −1 −1 1 ] [
G H H G
]−1 ·  100 1  = 2(a+ 2b+ c)κ+ 2
4(ac− b2)κ2 + 2(a+ c)κ+ 1
(19)
Based on the definitions of {a, b, c, d}, by using Cauchy inequality and completing square, we could prove: ac− b2 ≥ 0, a+ c ≥ 0, a+ 2b+ c ≥ 0, and a− 2b+ c ≥ 0. Treat κ ∈ (0,∞) as a variable. The derivative of (19) with respect to κ is a fraction. Its denominator is always larger than zero, while the numerator is a quadratic polynomial:
8(b2 − ac)(a+ 2b+ c)︸ ︷︷ ︸ ≤0 κ2 + 16(b2 − ac)︸ ︷︷ ︸ ≤0 κ− 2 (a− 2b+ c)︸ ︷︷ ︸ ≥0
(20)
Thus ρ1 + ρ4 − ρ2 − ρ3 monotonically decreases when κ increases in (0,∞). Its minimum value equals 0 and is attained when κ→∞. Our proof completes.

heading: A.5 UNIMODAL HIGH-DIMENSIONAL CASE
text: Without loss of generality, we consider the case when WWT = I. We emphasize this is a mild assumption. Because when W doesn’t satisfy this condition, we could perform SVD decomposition W = UΣV and convert the FC layer Wh + b into three sequential linear layers: h → Vh → ΣVh→ UΣVh+b, where now the weight matrix in the last layer satisfies UUT = I. Under this assumption B = (WWT )−1 = I, (8) could be further simplified as:
Wx = [( 1
γ − 1)J′,2 + I]−1 = (κJ′,2 + I)−1 (21)
where J′ = J(0) is known, and for simplicity, we have denoted κ = 1/γ − 1. From this expression, it is clear that the second term of (5) guarantees numerical stability, since J′,2 is not invertible and without the identity matrix (21) doesn’t exist. Before moving on, we define a helper vector eij ∈ RK , whose i-th and j-th entries equal 1 and −1, respectively, and all other entries equal 0. Then considering the relation of z′ = Wxz, we have:
(z′i − z′j)(zi − zj) = (z′,Teij)(eTijz) = z′,Teije T ijW −1 x z ′ (22)
This implies that if eijeTijW −1 x is a positive semi-definite (PSD) matrix for arbitrary i 6= j, then the order preserving property holds (i.e., if zi > zj , then z′i > z ′ j). Since eije T ijW −1 x is asymmetrical, examining whether it is PSD requires us to focus on the summation of its transpose and itself. Namely, eijeTijW −1 x is PSD if and only if the eigenvalues of [eije T ijW −1 x + (eije T ijW −1 x )
T ] are all no less than zero. Notice that Wx and W−1x are symmetric as shown in (21), we have
eije T ijW −1 x + (eije T ijW −1 x ) T = eije T ijW −1 x + W −1 x eije T ij
= eije T ij(κJ ′,2 + I) + (κJ′,2 + I)eije T ij
= 2eije T ij + κ(eije T ijJ ′,2 + J′,2eije T ij)
(23)
For simplicity we denote the set of κ which can make all eigenvalues of the aforementioned matrix non-negative as: Γij = {κ > 0|λmin(2eijeTij + κ(eijeTijJ′,2 + J′,2eijeTij)) ≥ 0}, where λmin(·) represents the minimum eigenvalue of a matrix. Then a necessary and sufficient condition for the order persevering property is κ ∈ ∩i 6=jΓij . An immediate corollary is: if we want the final prediction to change after adding Wx, then we must have 1γ − 1 /∈ ∩i6=jΓij .

heading: B EXPERIMENT DETAILS AND ABLATION STUDY

heading: B.1 DEFINITION OF OMEGA
text: Take AV-MNIST as an example. We have defined the following types of perturbations.
• Gaussian noise: Perturb every element x in the image (or audio) feature by a Gaussian noise in the form of xnoise = (1 + )x, where ∼ N(0, ω20).
• Missing entries: For the audio modality, we randomly select ω1 consecutive columns (or rows) and all elements in there are set to 0. This corresponds to missing time frames (or frequency components).
• Bias noise: For the image modality, we randomly select an image patch with size of ω2 × ω2, and every element x in there is perturbed by a given amount: xnoise = (1 +ω3)x. This corresponds to change of illumination. (see Figure 5).
• FGSM attack: We use the Cleverhans libarary (Papernot et al., 2016) to implement a FGSM attack on the input image/audio feature with l∞ norm. The maximum overall norm variation is ω4 (e.g., ω4 = 0.3).
• PGD attack: Cleverhans is adopted to implement a PGD attack on the input image/audio feature with l∞ norm. The maximum overall norm variation also equals ω4, and that for each inner attack is ω5 (e.g., ω5 = 0.01), and the number of inner attacks is ω6.
Unless explicitly mentioned, ω6 is set to 20.

heading: B.2 EXTRA RESULTS ON AV-MNIST
text: Impact of hyper-parameters Following Figure 4, we plot model accuracy versus magnitude of noise/attack with different gamma values in Figure 6. For comparison purpose, Figure 4 is re-drawn in its first row. As demonstrated in Figure 6, alongside the increasing of γ, the gap between the orange and blue lines becomes smaller. However, the orange lines consistently appear upon the blue lines, indicating that our method works in a wide range of gamma value and that hyper-parameter tuning is relatively easy in our method. We emphasize that γ = 1.0 is equivalent to trivial fusion.
Perturbation on the image modality Unlike in the main text, here we consider perturbations on the image modality, while the audio modality is assumed clean. Consequently, we invoke our Jacobian regularization method for the image modality, while keep the audio modality unchanged. We have deliberately enlarged the magnitude of the noise/attack, and results are reported in Table 4.
As shown in Table 4, confidence-weighted summation still doesn’t outperform purely mean fusion in all cases. Regards to statistical fusion, the second column demonstrates that this time our Jacobian regularization might lead to accuracy drop on clean data. The phenomenon that robust network
is less accurate on clean data might happen as suggested by Tsipras et al. (2019). However, we notice that such phenomenon doesn’t occur in Table 1. We believe that such difference might lie intrinsically in modality difference, and we will explore it in our future work.
Noise on both modality We further consider the setting when both image and audio modality are perturbed by Gaussian noise in Table 5, which differs ours from Kim & Ghosh (2019).

heading: B.3 EXTRA RESULTS ON RAVDESS
text: Impact of hyper-parameters Following Figure 6, we also plot model accuracy versus magnitude of noise/attack with different gamma values in the emotion recognition experiment. The results are shown in Figure 7. The multimodal network with our Jacobian regularization invoked (i.e., the orange line) almost surpasses the network without that invoked (i.e., the blue line) in all cases.
Perturbation on the image modality Here we consider perturbations on the image modality. Similarly, we have deliberately enlarged the magnitude of the noise/attack, and results are presented in Table 6. As demonstrated in the table, the mutimodal network enhanced by our method usually achieves the best accuracy among all models except in few cases. Also, our Jacobian regularization could usually improve model accuracy except in the clean case. We notice that in some extreme cases when the image modality is severely perturbed, the multimdoal network might be even worse than the regular unimodal audio network. This, in turn, implies that multimodal fusion is not always the winner, especially when large noise/attack is presented.

heading: B.4 EXTRA RESULTS ON VGGSOUND
text: Impact of hyper-parameters We also plot model accuracy versus magnitude of noise on video modality with different gamma values in the VGGSound experiment in Figure 8.
Unknown Noise on Single Modality To capture another setting, here we consider that the video modality (or audio modality) is corrupted, and both modalities are trained with robust add-ons in Ta-
ble 7 and 8, respectively. We emphasize that this setting is valid because in real-world applications, we sometimes do not know which modality is corrupted. As shown in these tables, the multimodal networks fused by our method could usually achieve higher accuracies in most cases, except when
both unimodal networks are learned with free-m training 3. Moreover, the fact that our Jacobian regularized multimodal network with a corrupted modality still outperforms the unimodal network with clean data in all cases demonstrates the positive effect of an extra modality.
3We neglect the cases when networks are tested on clean data, since it is possible that a robust network is less accurate on clean data as suggested by Tsipras et al. (2019).

Review summaries:
This paper propses a late-fusion algorithm for multimodal learning. This algorithm serves to improve the robustness against adversarial attacks and random corruptions. Assuming that which modality is perturbed, this paper propses to leverage Jacobian regularization and conditional independence assumption to fuse predictions from different modalities. Moreover, this paper also provides rigorous error bounds on it error and explain the necessity of extra modality.  In summary, this paper's main contributions are as follows:  - Propose a late-fusion algorithm based on Jacobian regularization and conditional independence; - Derive theoretical error bounds and demonstrate the biasing effect of the extra modality; - Conduct comprehensive evaluation on adversarial attacks and random corruptions and outperform baseline late-fusion algorithms.

The paper proposed a training-free robust multimodal learning late-fusion methods via sample-wise Jacobian regularization. The key idea of the work is to minimize the Frobenius norm of a Jacobian matrix, so that the multimodal prediction is stabilized. The paper demonstrate the good efficacy on both adversarial attacks and random corruptions setting on multimoda datasets such as AV-MNIST/VGGSound/RAVDESS. 

In this paper, the authors proposed a training-free late-fusion method for robust multimodal learning. They specifically considering its performance under adversarial attacks and random corruptions which usually confuse the model by introducing noise to the input data. To promote the multimodal prediction robust to attacks, they propose to minimize the Frobenius norm of Jacobian matrix so that the prediction becomes stable to the perturbation of inputs. They also provide a theoretical error bound of their method. The experimental results outperfom other late-fusion methods.