{'fraction_non_alphanumeric': 0.043614825450955755, 'fraction_numerical': 0.013461365879924617, 'mean_word_length': 4.3823860891965865, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Training-Free Robust Multimodal Learning via Sample-Wise Jacobian Regularization

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 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.


Introduction
Deep fusion models have recently drawn great attention of researchers in the context of multimodal learning [30,1,24,32,33] as it provides an easy way to boost 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 [15]. Over the past few years, numerous multimodal fusion methods have been proposed at different levels: early-, middle-, and late-fusion [6]. 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 [26]. 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 late-fusion 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 training data) and thus are convenient to deploy.
Although late-fusion is a mature topic in the literature, its performance under adversarial attacks [20,27] and random corruptions [37,15] is rather underexplored. In this paper, we address the problem of robust late-fusion by utilizing Jacobian regularization [29,12,11,3] and conditional independence assumption [26]. 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 [13], and the solution is obtained with only a little computational overhead.  Fig. 1. Illustration of the proposed robust late-fusion method. Before unimodal raw logit zA and zB are fused, the Jacobian regularization technique is applied. Roughly speaking, it will enforce the derivative of p with respect to zA becomes smaller. Thus, when zA is perturbed to z A (due to random corruption or adversarial attack), the change of multimodal prediction ||p −p|| will be limited to some extent. For illustration purpose, all variables here (e.g., p, zA) are drawn in one-dimensional space.

Preliminary
Network Robustness To verify 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 [9,20,2] and (ii) random corruptions such as Gaussian noise, missing entries or illumination change [37,15]. Correspondingly, many methods have been proposed to offset the negative effects of such perturbations. Adversarial training based on projected gradient descent [20] is one strong mechanism to defend against adversarial attacks, and the recent Free-m method [25] is proposed as its fast variant. Besides adversarial training, several regularization techniques are also proven to have such capability, such as Mixup [36], Jacobian regularization [12]. Alternatively, with regard to random corruptions, Mixup and Jacobian regularization are also effective in this case [35,11]. Another powerful approach is stability training [37], 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 [6,15], semantic segmentation [5,7], audio recognition [8,4]. Based on where the information is exchanged between different modalities, multimodal fusion methods could be classified into three kinds: (i) early-fusion [31,6], (ii) middle-fusion [15,32], and (iii) late-fusion [6,18]. 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, [22], [28], and [14] exploited gating networks to deal with random corruptions, adverse or changing environments. Afterwards, [15] 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.

Sample-Wise Jacobian Regularization
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 represents input features and y ∈ {0, 1, · · · , K − 1} represents the true label. We train two unimodal networks separately, each corresponding to one modality. Given a specific input feature x A , the first unimodal network calculates the class prediction p A ∈ R K by: where σ 1 (·) represents the Softmax function, z A ∈ R K represents the raw logit, and h A = f A (x A ) ∈ R H represents the feature being fed into the last layer 6 .
Here W A ∈ R K×H and b A ∈ R K are the learnable weight and bias of the last linear layer, respectively. Similarly, the second unimodal network provides the class prediction Based on the conditional independence assumption [16], the basic statistical late-fusion method generates the final class prediction as [6]: 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 ∈ R K contains the occurring frequencies of each class, calculated from the training dataset.
Our proposed approach builds upon (2). Specifically, we consider adding two K × K weight matrices {W a , W b } ahead of {z A , z B }, respectively, right before they get activated. Consequently, the final multimodal prediction is re-calibrated as: 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 W b , 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 W a , we first calculate the derivative of p with respect to h A (See Supplementary): Then we minimize the following regularized Jacobian loss with respect to W a : 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 h A 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 {x A , x B }, once W a is determined, so are p and J via (3) and (4), respectively. Thus, (5) is well-determined and non-linear with respect to W a .
We propose a heuristic iterative method making the above optimization problem tractable in Algorithm 1. Our key is to decouple J from W a . Namely, in step 5 of Algorithm 1, all terms are known except W a , and thus the relaxed loss is convex. After writing out all terms of ∂L (t) /∂W a , we observe that it is a Sylvester equation [13]. It has a unique solution and the run time is as large as inverting a K × K matrix and hence affordable. See Supplementary 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 t max = 1 (i.e., p = p (1) ). Second, if we know data from modality B are also perturbed, we could solve W b in a similar manner. Third, 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 Supplementary for details). Some immediate corollary include: (i) when ∼ N (0, Σ), the bound is simplified to E[||p ,noise − p ||] ≤ l ( γK 2(1−γ) ) 1/2 Tr[Σ], (ii) when the L 2 norm of is constrained smaller than δ (usually assumed in adversarial attacks), the bound is simplified to ||p ,noise − p || ≤ l δ ( γK 2(1−γ) ) 1/2 .

Necessity of the Extra Modality: Biasing Effect
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 W x before the raw logit z being fed into the Softmax activation, re-calibrating the final unimodal prediction as: Here W x can be solved by analogy to (5). However, we observe that the introduction of W x usually won't change the final prediction. For an intuitive understanding, we consider the TwoMoon example used by [33]. 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.
(a) Train on clean data

97.00%
------ 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 W x 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 W x in unimodal binary classification. See Supplementary for rigorous mathematical proof. Moreover, in a K-class (K > 2) classification problem, the final prediction might change if γ is sufficiently small. 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 γ ↓).
(a) Train on clean data  Lemma 1. Statistical fusion is equivalent to concatenating the raw logits: 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 Supplementary). Now assume that the sample-wise W a is added, then the final multimodal prediction becomes σ 1 (W a z A + z B − ln freq). Alternatively, in the unimodal context, our method adds a sample-wise W x and the final prediction becomes σ 1 (W x z). Comparing these two expressions, we observe that the introduced W a and W x occur ahead of the perturbed modality. In the unimodal case, the Softmax function σ 1 directly applies onto W x z, while in the multimodal case, the extra modality information z B acts as a bias term inside the Softmax. Since entries of z B have different values, it might impact the final prediction.

AV-MNIST
AV-MNIST [24,30] 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 [24]. 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 [30], we use LeNet5 [17] and a 6layer 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 Supplementary 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 [20], (iii) free-m training (freT) by [25], (iv) stability training (stabT) by [37], and (v) Mixup (mixT) by [35]. 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 [34] and adapt 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 [15], our experiments are mainly conducted on the case when one modality is perturbed. Results on various other settings are reported in Supplementary.
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 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. The first row of 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. We also plot model accuracy versus magnitude of noise/attack with different gamma values in the second and third row of Figure 4. 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.
On the other hand, 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 2. As shown in Table 2, 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 [27]. 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.

Emotion Recognition on RAVDESS
We consider emotion recognition on RAVDESS [19,33]. 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 3.
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 (e.g., perturbation on the image modality) are presented in Supplementary.  We consider classification task on a real-world data set VGGSound [4], 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 audiovideo 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 [10] 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 4. 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. In Supplementary, we capture another important setting: the video modality (or audio modality) is corrupted, but we don't know which one. Thus, both modalities are trained with robust add-ons. Please refer to Supplementary for more details.

Discussion and Future Work
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 W x , we deliberately add another bias term b x 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. In this subsection, we first prove Lemma 1 shown in the main text. For conciseness, we denote the k-th element of p A as p A,k . Similar notations are adopted for p B , z A , z B , and ln freq. The k-th element of the left-hand side (LHS) in Lemma 1 can be simplified as: 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 (S1), we use the definition p A = σ 1 (z A ), p B = σ 1 (z B ), and have introduced two normalization constants T A = || exp(z A )|| 1 and T B = || exp(z B )|| 1 . As a result, an equivalent form of (3) is p = σ 1 (W a z A + W b z B − ln freq). Based on the chain rule and z A = W A h A + b A , we have: where we have used the Jacobian matrix of Softmax function to calculate J . We emphasize that J is symmetric.

A.2 Solution of Regularized Jacobian Loss
Because Frobenius norm is convex and a weighted summation with positive coefficients preserves convexity, the relaxed optimization problem shown in Step 5 of Algorithm 1 is a convex function with respect to W a . This implies that any local minimum is also the global minimum. We calculate the derivative of the relaxed optimization problem with respect to W a and set it to zero, yielding: where for simplicity we have defined A = (1/γ − 1)J (t),2 ∈ R K×K and B = (W A W T A ) −1 ∈ R K×K . (S3) is known as Sylvester equation 7 , and it has been proven that if A and −B do not share any eigenvalues, then the Sylvester equation has an unique solution. 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 wellgeneralized neural network normally doesn't have singular weight matrices [21]. Thus, 0 is almost surely not the overlapping eigenvalue and a unique solution of (S3) is guaranteed. As demonstrated by [13], the matrix equation (S3) 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 (5) guarantees numerical stability.

A.3 Proof of Theorem 1
In the last iteration t = t max → ∞, noticing W a satisfies (S3), we have: and: Thus we have: (S6) where we have taken advantage of the fact that W a W T a is a positive semi-definite matrix, and thus its trace is no less than zero. Based on the first-order Taylor expansion, we have: Thus we have: where we have utilized the l-Lipschitz continuity of f A .

A.4 Unimodal 2D Case
Consider an illustrating example when K = 2 (a.k.a. unimodal binary classification). In this case, the final prediction is two dimensional: p = [p 1 p 2 ] T . We denote the raw logit after adding the W x matrix as z = W x z, and its entries are denoted as z 1 and z 2 . Let us begin by explicitly writing out the first line of (S3). After some simplification, we obtain: where κ = 1/γ−1. Here we have explicitly written out the elements of WW T and merged some terms of J ,2 into a = p 1 p 2 (w 2 11 + w 2 12 ), b = p 1 p 2 (w 21 w 11 + w 12 w 22 ), and c = p 1 p 2 (w 2 21 + w 2 22 ). The unknowns occur in W x : After some algebra, we can solve all unknowns via the following block matrix inversion:  where for simplicity, we have defined: Noticing that z = W x z, we can prove that: 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 z 1 > z 2 , then z 1 > z 2 and vice versa. This order preserving property is equivalent to say that adding the W x 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: Based on the definitions of {a, b, c, d}, by using Cauchy inequality and completing square, we could prove: ac − b 2 ≥ 0, a + c ≥ 0, a + 2b + c ≥ 0, and a − 2b + c ≥ 0.

A.5 Unimodal High-Dimensional Case
Without loss of generality, we consider the case when WW T = 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 UU T = I. Under this assumption B = (WW T ) −1 = I, (S3) could be further simplified as: 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 (S16) doesn't exist. Before moving on, we define a helper vector e ij ∈ R K , whose i-th and j-th entries equal 1 and −1, respectively, and all other entries equal 0. Then considering the relation of z = W x z, we have: This implies that if e ij e T ij W −1 x is a positive semi-definite (PSD) matrix for arbitrary i = j, then the order preserving property holds (i.e., if z i > z j , then z i > z j ). Since e ij e T ij W −1 x is asymmetrical, examining whether it is PSD requires us to focus on the summation of its transpose and itself. Namely, e ij e T ij W −1 x is PSD if and only if the eigenvalues of [e ij e T ij W −1 x + (e ij e T ij W −1 x ) T ] are all no less than zero. Notice that W x and W −1 x are symmetric as shown in (S16), we have For simplicity we denote the set of κ which can make all eigenvalues of the aforementioned matrix non-negative as: Γ ij = {κ > 0|λ min (2e ij e T ij + κ(e ij e T ij J ,2 + J ,2 e ij e T ij )) ≥ 0}, where λ min (·) represents the minimum eigenvalue of a matrix. Then a necessary and sufficient condition for the order persevering property is κ ∈ ∩ i =j Γ ij . An immediate corollary is: if we want the final prediction to change after adding W x , then we must have 1

B.1 Definition of Omega
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 x noise = (1 + )x, where ∼ N (0, ω 2 0 ). -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: x. This corresponds to change of illumination. (see Figure 5). -FGSM attack: We use the Cleverhans libarary [23] 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.

B.2 Extra Results on RAVDESS
Impact of hyper-parameters Following Figure 4 of the main text 8 , 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  6. 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 5. 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.

B.3 Extra Results on VGGSound
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 7.
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 Table 6 and 7, 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 9 . 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.