Few-shot Multi-Task Learning of Linear Invariant Features with Meta Subspace Pursuit
Chaozhi Zhang Lin Liu Xiaoqun Zhang


Hong-Ou-Mandel Interference in a temporal-average-inversion-symmetric chain
    Zhoutao Lei
    September 9, 2024
===========================================================================



§ ABSTRACT 
Data scarcity poses a serious threat to modern machine learning and artificial intelligence, as their practical success typically relies on the availability of big datasets. One effective strategy to mitigate the issue of insufficient data is to first harness information from other data sources possessing certain similarities in the study design stage, and then employ the multi-task or meta learning framework in the analysis stage. In this paper, we focus on multi-task (or multi-source) linear models whose coefficients across tasks share an invariant low-rank component, a popular structural assumption considered in the recent multi-task or meta learning literature. Under this assumption, we propose a new algorithm, called  (abbreviated as ), that provably learns this invariant subspace shared by different tasks. Under this stylized setup for multi-task or meta learning, we establish both the algorithmic and statistical guarantees of the proposed method. Extensive numerical experiments are conducted, comparing  against several competing methods, including popular, off-the-shelf model-agnostic meta learning algorithms such as . These experiments demonstrate that  achieves superior performance over the competing methods in various aspects.




§ INTRODUCTION



Modern machine learning methods generally demand extensive amount of data to achieve good learning performance. However, in certain application domains such as clinical medicine and social sciences, collecting a large amount of training data can be extremely cost and time consuming, and sometimes even unethical due to privacy concerns. As a result, data scientists and statisticians will have to face the so-called “data-scarce” or “data-hungry” challenge when deploying those modern machine learning methods <cit.>. 

To bypass this difficulty, the multi-task or meta learning paradigm <cit.> has emerged as a promising option both in theory and in practice <cit.>, by borrowing information across different data sources/tasks but sharing certain invariances. Intuitively, multi-task or meta learning methods are advantageous because they are designed to acquire those invariant information shared by all tasks during the learning process. 



To make progress in our theoretical understanding of multi-task or meta learning, a large body of the current literature focuses on the multi-task large-dimensional linear models <cit.>, a stylized statistical model but still preserving some essential features of more general multi-task or meta learning problems encountered in real life. Specifically, we are given a set of T tasks, indexed by [T] {1, ⋯, T}. For each task t ∈ [T], we are given access to m independent and identically distributed (i.i.d.) observation pairs (_t∈^m × d, _t∈^m), where _t and _t, respectively, denote the design matrix of d features and the scalar-valued responses/labels, over all m samples in task t. The total sample size is denoted as N ≡ T m. We further assume that for every task t, the response _t follows a noise-corrupted linear model:

    _t = _t _t^∗ + _t, t = 1, ⋯, T,

where _t ∈ℝ^m represents the exogenous white noise associated with task t. y_t, j, _t, j, _t, j are reserved for the response, the features, and the random noise for sample j in task t. In addition, x_t, j, k is the k-th coordinate of the feature vector _t, j. We let Θ^∗ [_1^∗, ⋯, _T^∗]^⊤∈^T × d be the collection of the true regression coefficients over all tasks. Since it becomes more difficult to estimate Θ when d becomes larger relative to m (corresponding to data-scarcity/data-hungry regime <cit.>), further structural assumptions on the model (<ref>) are indispensable. A popular choice is the so-called “Hard Parameter Sharing (HPS)” condition by assuming that the regression coefficients over all tasks enjoy a factor structure Θ^∗^∗^∗, where ^∗∈^T × s is the task-varying low-rank component with rank s ≤ d and ^∗∈^s × d is the task-invariant component that maps the intrinsic task-varying coefficient _t^∗∈^s linearly to _t^∗∈^d, for every t ∈ [T]. Intuitively, once the task-invariant ^∗ is at our disposal, the original high-dimensional task can be reduced to easier, lower-dimensional task by projecting the features onto {^∗}, the linear span of the column space of ^∗. 



 §.§ Related Work


Although exceedingly simple, the model (<ref>) has been analyzed extensively as a first-step towards understanding multi-task or meta learning. Due to space limitation, we only highlight several works that are relevant to our development. <cit.> introduced a Method-of-Moments () algorithm that learns the task-invariant subspace {^∗} via Singular-Value-Decomposition (SVD) of a y^2-weighted empirical Gram matrices of the features , where we recall that y is the response variable. The  method, although extremely simple, is nearly rate optimal for estimating {^∗} in sine-angle distance by matching the minimax lower bound derived in the same paper <cit.>, which scales roughly in order √(d s/T m), where we recall that d, s, m, T are respectively the feature dimension, the rank of the task-invariant subspace ^∗, the sample size of each individual task, and the total number of tasks. Lower bounds in a similar spirit have also been derived in <cit.> and in <cit.>, where the authors studied a variant of the HPS linear model by allowing the low-rank subspace/representation ^∗ to slightly vary across tasks. 

Despite being rate-optimal, as shown later in Section <ref>,  does not perform well in practice, even in very standard settings where the feature is drawn from isotropic Gaussian. A modified version of  in <cit.> offers improved performance in estimating the subspace, particularly when data is extremely scarce. However, this modified  struggles to perform that well in practice as the feature dimension and sample size grow (see Figure <ref>). Another alternative method presented in <cit.> incorporates the Burer-Monteiro factorization method <cit.>. However, this method fails to effectively handle cases with very scarce data, and the derived rate of estimating ^∗ does not shrink to zero as the number T of tasks grows, hence a suboptimal rate compared against the minimax lower bound √(d s/T m). <cit.> proposes an algorithm that involves alternating minimization and gradient descent techniques. Nevertheless, as shown in Figure <ref>, this method still requires a substantial amount of data, which may be difficult to acquire in some applications.

In a slightly different vein, <cit.> demonstrate that Model-Agnostic Meta-Learning (MAML) <cit.> and Almost-No-Inner-Loop algorithm (ANIL) <cit.>, which are not designed for linear models, are still capable of learning the shared representation, or the task-invariant subspace ^∗. However, the outcomes of these algorithms are not well-suited for data scarcity situations, as their analysis relies on drawing independent mini-batches at each iteration. As shown in our numerical experiments, ANIL requires more data to achieve similar error level(see Figure <ref>). These empirical performance suggests that when the data generating mechanism satisfies certain complexity-reducing structures (e.g. the HPS linear model), MAML and ANIL could be sub-optimal because they fail to leverage such structures at the expense of being model-agnostic.

Finally, we mention that <cit.> proposed a trace-norm regularized estimator and derived error bounds of their proposed method[Their error bounds were later improved to nearly optimal by <cit.> under Gaussianity with Gordon's Convex Gaussian Min-Max Theorem [CGMT] <cit.>. However, <cit.> further imposes the constraint T < m; see the statement of their Theorem 3.]. As will be evident later, this approach is the closest to ours as the trace norm can be viewed as a relaxation of the matrix rank used in our method. Nevertheless, they do not provide a satisfactory algorithmic solution to address the optimization problem[In their numerical experiments, the BFGS method is directly used to solve the relevant optimization problems, which may lead to sub-optimal solution. In Section <ref>, we elaborate on this issue when evaluating the finite-sample performance of different methods.].



 §.§ A Prelude to Our Approach and Our Contributions

Overall, the existing methods present various strengths and weaknesses, making it essential to develop new approaches that can effectively handle data scarcity by learning the task-invariant subspace effectively. In this paper, we propose a method of directly learning the shared representation, or equivalently the invariant subspace in the HPS linear model, inspired by the matrix rank minimization problem. We will formally introduce our new algorithm later in Section <ref>. In this section, we provide some intuitive explanation of our proposed algorithm.

Specifically, we address the problem of learning task-invariant subspace by treating it as a matrix rank minimization problem. The corresponding constrained matrix rank minimization problem is formulated as follows:

    min_X     rank (X) 
    s.t.     𝒜 (X) = b,

where X is the unknown matrix, 𝒜 is a linear map and b is a vector. As in <cit.>, there is a more robust formulation of (<ref>):

    min_X     ‖𝒜 (X)-b ‖_2 
    s.t.     rank (X) ≤ s,

which we will really consider in our methods.

The matrix rank minimization problem is closely related to the compressed sensing problem. They become equivalent when X is reduced to an unknown vector, and the objective function is replaced by ‖ X ‖_0. Solving compressed sensing and matrix rank minimization problems exactly, in the form (<ref>), are well known NP-hard problems <cit.>. To find an approximate solution to the compressed sensing problem, it is common to replace the objective function with the ℓ_1 norm, which is the convex envelope of the ℓ_0 norm. Various algorithms, such as iterative hard thresholding <cit.>, subspace pursuit <cit.>, compressive sampling matching pursuit <cit.>, iterative thresholding with inversion <cit.>, and hard thresholding pursuit <cit.>, are typically employed for the original ℓ_0 minimization problem. For the ℓ_1 optimization problem, several algorithms, such as those introduced in <cit.>, <cit.>, <cit.>, and <cit.>, are commonly used.

In analogy, nuclear norm minimization, which replaces the rank of the matrix with its nuclear norm, solves the matrix rank minimization problem (<ref>) approximately. The method proposed in <cit.> exemplifies this approach. As the nuclear norm minimization problem is equivalent to some semi-definite programming (SDP) problem, interior-point methods can be used to solve it, as shown in <cit.>. Several other first-order algorithms, including singular value thresholding (SVT), fixed point continuation (FPC), and the approximate SVD-based FPC algorithm (FPCA), have also been proposed in <cit.> and <cit.>. In our paper, we opt to utilize the iterative hard thresholding algorithm mentioned in <cit.>. This algorithm is a generalization of the iterative hard thresholding algorithm used in compressed sensing <cit.>. Within the context of meta-learning, we designate this algorithm as  (). We substantiate its efficacy by demonstrating that it can attain a convergence outcome analogous to the one reported in <cit.>. In practical applications, our approach demonstrates improved statistical properties and computational efficiency empirically.

We now summarize our main contributions as follows:


  * On the methodological side, we develop a new, iterative algorithm  for learning the task-invariant subspace in the multi-task HPS linear models, contributing a new method to the multi-task/meta/invariant learning paradigm <cit.>.  is easy and transparent to implement, .

  * On the theoretical side, we establish how fast the regression coefficients  and task-invariant subspace  outputted by  converge to the truth ^∗ and ^∗, as m, d, s vary. In particular, we directly analyze the iterations of  and establish the convergence rates at k-th iteration, for every k ≥ 1. As a consequence, our analysis directly reveals the iteration complexity of .

  * Empirically, through extensive experiments based on simulated and real datasets, we demonstrate that  outperforms most of the competing methods on several aspects. Several future directions for theoretical investigation are also hinted by the results of our numerical experiments. For more details, see Section <ref>. 




 §.§ Organization of the Paper

The rest of our paper is organized as follows. In Section <ref>, we introduce the basic problem setup, the statistical model being analyzed, and key underlying assumptions. The new  algorithm is then proposed in Section <ref>, with its theoretical properties, in terms of sample and computational complexities, established in Section <ref>. In Section <ref>, we conduct numerical experiments with both simulated and real datasets that demonstrate the performance of our proposed method in practice, together with an extensive comparison between our method and other competing methods. Finally, we conclude our paper in Section <ref> with a discussion on future research directions.



 §.§ Notation


Before proceeding, we collect some notation frequently used throughout the paper. The distributions for multivariate Gaussian and sub-Gaussian random variables are denoted as (μ,Σ) and (μ,Σ), respectively, where μ is the mean and Σ is the covariance matrix. For sub-Gaussian random variables , by (μ, Σ), we mean that the population covariance matrix [( - μ) ( - μ)^⊤] is dominated by Σ in the positive semidefinite sense. In general, we reserve bold fonts for vectors and matrices and regular fonts for their elements.

Given a real-valued vector 𝐛, let 𝐛_q denote its ℓ_q norms. For a real-valued matrix , _2, _, _∗, and [] represent its spectral norm, Frobenius norm, trace norm, and rank respectively. () and ^⊥ () denote the column space spanned by  and its orthocomplement. With slight abuse of notation, we denote ^⊥ as a matrix with (^⊥) = ^⊥ (). In general, ^⊥ is not unique, and we just pick one arbitrarily. Throughout the paper, we use the sine angle distance to measure the distance between the column spaces spanned by two different matrices _1 and _2, defined as sin∠ (_1, _2)‖_1^⊤_2^⊥‖_2 <cit.>. Given any positive integer k, we let _k denote the k × k identity matrix. Given a matrix ∈^k_1 × k_2, we introduce the following vectorization operator to flatten the matrix :

    (⃗) [ 𝐦_1, ·^⊤, ⋯, 𝐦_k_1, ·^⊤]^⊤.

If  is a symmetric positive semidefinite matrix, let ^† denote its Moore-Penrose pseudo inverse (also known as the ridgeless regularized inverse ^†lim_λ↓ 0 ( + λ)^-1): let = ^⊤ be the eigen-decomposition of  and let  be the diagonal matrix of only the non-zero eigenvalues and  be the corresponding eigenvectors, then ^†≡^-1. Finally, given two matrices  and  of appropriate sizes, ⊙ denotes their matrix Hadamard product.



§ PROBLEM SETUP AND MAIN ASSUMPTIONS


    





 §.§ The statistical model



As alluded to in the Introduction, the Hard Parameter Sharing (HPS) linear model posits a factor structure of the regression coefficients among different tasks: There exists a task-invariant matrix ^∗∈^s × d and a task-varying, but lower-dimensional matrix ^∗∈^T × s, with s ≤ d, such that Θ^∗ = ^∗^∗. Here ^∗≡ [_1^∗, ⋯, _T^∗]^⊤ can be interpreted as the intrinsic task-specific coefficients and ^∗ represents the task-invariant linear map that lifts _i∈^s to _i∈^d for all i ∈ [T]. ^∗ is sometimes called the task-invariant linear representation in the literature <cit.>. ^∗ is assumed to have orthonormal columns as we are only interested in the linear space {^∗} spanned by the columns of ^∗. To avoid clutter, ^∗ will also mean {^∗} when it is clear from the context.

In multi-task or meta learning, an essential idea is to learn the low-dimensional shared representation ^∗ by  trained from multiple data sources. Based on learnt , the original high-dimensional problem is reduced to lower dimensions, thus enjoying improved sample and computational complexity. Following the current literature <cit.>, the quality of learning ^∗ by  is gauged by the sine angle distance sin∠(, ^∗).

With Model (<ref>) in place, we also need to specify some assumptions on the data. First, we impose the following distributional assumptions on the covariates _i and the additive noise _i:

For every task t ∈ [T], x_t, j, ki.i.d.∼ (0, 1) for j = 1, ⋯, m and k = 1, ⋯, d; and ε_t, j∼ (0, σ^2) for j = 1, ⋯, m.


The distributional assumption imposed on _t implies that it satisfies the Restricted Isometry Property (RIP) (or properties alike such as the Restricted Eigenvalue Conditions <cit.>) with high probability, a common condition imposed on the covariates to ensure identifiability (i.e. well-posedness along any sparse directions). Similar to the existing literature, we also assume the following condition on the intrinsic task-specific regression coefficients ^∗:
[Task Diversity]
Let λ_1, ⋯, λ_s be the largest to the smallest eigenvalues of the Gram matrix of task-varying components averaged over tasks, Ξ^∗ T^-1^∗⊤^∗≡ T^-1∑_t = 1^T_t^∗_t^∗⊤. There exists some universal constant L > 0 such that λ_s≥ L.

This assumption essentially requires that the task-varying information in the regression coefficients is sufficiently diverse across tasks, to the extent that their Gram matrix Ξ^∗ averaged across tasks is well-conditioned.

To ease exposition, we introduce some further notations: 𝕏∈^mT × dT denotes the block-diagonal matrix with block-diagonal elements _t∈^m × d for t = 1, ⋯, T, = (_1^⊤, ⋯, _T^⊤)^⊤∈^N × 1 and = (_1^⊤, ⋯, _T^⊤)^⊤∈^N × 1, with N  m T. Additionally, given any Θ∈^T × d, we define the linear mapping : ^T × d→^N × 1 as follows

    (Θ) 𝕏·(⃗Θ) ≡[ [ _1      ;     ⋱   ;       _T;    ]]
        [ [ _1;  ⋮; _T;    ]].

Hence we can rewrite Model (<ref>) compactly in matrix form as:

    =  (Θ^∗) + .

Our formulation can also be adapted to the case where different tasks have different sample sizes m_1, ⋯, m_T. But to not lead readers astray, we decide to focus on the simpler, equal-sized case, i.e. m_1≡⋯≡ m_T≡ m.



 §.§ Restricted Isometry Property


Similar to the compressed sensing literature, In the pursuit of solving the aforementioned problem through matrix rank minimization, the concept of the Restricted Isometry Property (RIP) plays a pivotal role. While RIP for sparse vectors in compressed sensing was originally introduced in <cit.>, it has been extended to matrices to address matrix rank minimization, as seen in <cit.>.  Here, we define RIP in the context of Model (<ref>):

For any integer r where 1≤ r ≤ d, consider the linear operator :ℝ^T × d→ℝ^N × 1. The operator  satisfies the Restricted Isometry Property with the restricted isometry
constant δ_r(), where δ_r() is the smallest constant for which the following inequality holds:

    (1-δ_r())Θ_^2 ≤(Θ)_2^2 ≤ (1+δ_r())Θ_^2

for all Θ∈ℝ^T × d with rank(Θ)≤ r.


In our context, we define the linear operator =𝒜/√(m), where 𝒜 is as defined in Equation (<ref>). Equation (<ref>) can be rewritten as:

    (1-δ_r())(∑_t=1^T_t_2^2) ≤∑_t=1^T1/√(m)_t_t_2^2 ≤ (1+δ_r())(∑_t=1^T_t_2^2)

Since rank(Θ)≤ r, all _t lie within an r-dimensional subspace. Therefore, there exist an orthonormal matrix ∈ℝ^r × d(^⊤=_r) and a set of vectors {_t}_t=1^T such that _t=^⊤_t. Equation (<ref>) can then be equivalently expressed as:

    (1-δ_r())(∑_t=1^T_t^⊤_t)
        
        ≤∑_t=1^T1/m_t^⊤_t^⊤_t^⊤_t ≤
        (1+δ_r())(∑_t=1^T_t^⊤_t).

This version of the RIP condition can often be justified by the distributional assumptions on the feature vector , using matrix concentration inequalities <cit.>.


Let _1, ⋯, _m ∈^d_1 × d_2 be independent, centered random matrices of the same size, where _j = 0 and _j≤ L for every j = 1, …, m. Let ∑_j=1^m_j. Define the quantity ν(), which can be viewed as the matrix extension of the variance for scalar Bernstein inequality <cit.>, as follows:

    ν()    max{ (^⊤), (^⊤)}
       ≡max{ (∑_j=1^m_j_j^⊤), (∑_j=1^m_j^⊤_j)}.

Then we have:

    ℙ (≥ u) ≤ (d_1 + d_2) ·exp{-u^2 / 2/ν () + L u / 3} for all  u ≥ 0.



The proof of Lemma <ref> can be found in <cit.>. Armed with this lemma, we obtain the following:

The linear operator =𝒜/√(m) where 𝒜 is defined in Equation (<ref>), satisfies the Restricted Isometry Property with the restricted isometry constant:

    δ_r()≤√(8(a+1)r-4/3mlog2r/ϵ)

with probability at least (1-ϵ)^T(1-15/a^2)^rT for any a>0.


The proof of Theorem <ref> is presented in Appendix <ref>. It demonstrates that when the data generating mechanism follows Assumption <ref>, the associated linear operator 𝒜 satisfies the RIP condition with high probability. This guarantee forms the building block for the convergence analysis of our proposed method.



§ OUR METHOD



In this section, we are ready to describe the new methodology for learning the invariant subspace ^∗ under the HPS linear model. The method is coined as the  Algorithm, or  for short.


In a nutshell,  consists of two steps. In the first step, task-specific regressions are solved jointly to obtain the coefficient matrix estimator Θ̂. To simplify exposition, we assume that the rank s of Θ^∗ (or equivalently, the dimension s of the subspace ) is known, whence Θ̂ is the solution to the following constrained optimization problem:

    min_Θ     ℒ(Θ):=𝒜()-_F^2,
    s.t.     rank() ≤ s

Here, 𝒜 (Θ) is the linear mapping defined in Equation (<ref>). 
In the second step, with Θ just attained from step one, it is then straightforward to apply either QR matrix factorization or Singular Value Decomposition (SVD) to obtain the task-invariant subspace B. In the following, we introduce the proposed  () algorithm for solving this constrained optimization problem. 



 §.§ The  () Algorithm



The  algorithm can be viewed as a variant of the iterative hard thresholding (IHT) algorithm, initially introduced for solving compressed sensing problems in <cit.>. The IHT algorithm has since then been adapted and generalized to handle matrix rank minimization problems, similar to the formulation in Equation (<ref>). Here we further adapt the matrix rank minimization problem to solve the task-invariant subspace learning problem. Since (<ref>) and (<ref>) have the same form when 𝒜() and  are vectorized, we can use this kind of algorithm to solve (<ref>).

The  algorithm constitutes the following iterative subroutines:

    
  * Gradient Descent (GD) Step. In each iteration, a GD update is performed independently for each task. This step aims to reduce the task-specific loss. For the t-th task, the (k + 1)-th iteration is updated according to the following expression:
    
    _t^(k+1) = _t^(k) - ∇__tℒ () |__t = _t^(k) = _t^(k) + γ/m_t^⊤ (_t - _t _t^(k))

    where γ is the step size and _t^(k) is the regression coefficients updated after the k-th step to be defined later. The updated regression coefficients are also concatenated as Θ̂^(k+1) [_1^(k+1), ⋯, _T^(k+1)]^⊤.
    
    
  * Hard Thresholding (HT) Step. After the GD update, a hard thresholding operation is then applied to the singular values of the concatenated coefficient matrix Θ̂^(k + 1) just obtained. Specifically, we first use SVD to factorize Θ̂^(k + 1) as Θ̂^(k+1)=^(k+1)^(k+1)^(k+1)^⊤ where ^(k+1)∈ℝ^T × T, ^(k+1)∈ℝ^T × d and ^(k+1)∈ℝ^d × d are the matrices for the left singular vectors, the singular values, and the right singular vectors, respectively. Then, only the largest s singular values and the corresponding left and right singular vectors are retained, recalling that s is the task-invariant subspace dimension. This step enforces the rank constraint over iterative updates of . The result of the hard thresholding step is a new matrix Θ^(k + 1) = ^(k + 1)^(k + 1)^(k + 1)^⊤, where ^(k+1)∈ℝ^T × s and ^(k + 1)∈ℝ^d × s are the first s columns of ^(k+1) and ^(k + 1), and ^(k + 1)∈ℝ^s × s is the diagonal matrix with the first s diagonal singular values of ^(k + 1). For short, we denote this hard thresholding operation as _s (·), so _s (Θ̂^(k + 1)) = Θ^(k + 1).
    
    
  * Representation Learning Update. The right singular vectors ^(k + 1) ⊤ from the HT step is taken as the learnt task-invariant subspace ^(k + 1) at the (k+1)-th iteration, since it consists of orthonormal columns by construction.

The  Algorithm iteratively performs the above three steps until a stopping criterion is met. The complete procedure of the  is presented in Algorithm <ref>. It is noteworthy that the step size γ is an important hyperparameter – the algorithm converges to the desired solution only if γ is set appropriately.





 §.§ Theoretical analysis of 



In this section, our goal is to establish error bounds for the estimator Θ^(k) obtained by the  algorithm and the  ground truth Θ^∗. Additionally, we aim to constrain the sine angle distance between the subspaces 𝐁^(k) and ^∗.


In this context, we operate under the assumption that the Restricted Isometry Property (RIP) is valid, as established in Theorem <ref>. This assumption provides us the leverage to establish error bounds for the Θ^(k) estimator. Our primary theorem is presented below.


Under Assumptions <ref> to <ref>, we suppose the restricted isometry constant of =𝒜/√(m) is δ_s()(𝒜 is defined in Equation (<ref>) and δ_s() satisfies Theorem <ref>). If 1-1/2√(2)/1-δ_3s<γ≤1, for any a > 0, with probability at least (1 - ϵ) (1 - 15/a^2)^3s( 1 - 2/d - 14/m - 208/dm)^T, the following holds:

    Θ^∗ - Θ^(k)_F ≤[2√(2)(1-γ+γδ_3s())]^kΘ^∗ - Θ^(0)_F + 2√(2)/2√(2)(1-γ+γδ_3s())√(dTσ^2/m).

In particular, if γ=1,

    Θ^∗ - Θ^(k)_F ≤( 192(a+1)s-32/3mlog6s/ϵ)^k/2Θ^∗ - Θ^(0)_F + 2 √(2)σ/1 - √(192 (a+1) s - 32/3mlog6s/ϵ)√(dT/m).


Moreover, as k →∞, the above bounds can be simplified to:

    Θ^∗-Θ^(k)_F ≤ O(√(σ^2dT/m))


The proof of Theorem <ref> is deferred to Appendix <ref>. Our proof loosely follows the strategy of <cit.>. However, during the proof, we need to derive more refined error bounds depending on both the sample size m and the data dimension d. 

The bound in Frobenius norm shown in Theorem <ref> can then be used to establish convergence rate of the learnt task-invariant subspace ^(k) after k-th iteration of  to the truth ^∗.

We give a theoretically provable bound for γ according to Theorem <ref>. In practice, this bound could be wider.


Under the same assumptions of Theorem <ref>, we have:

    sin∠ (^(k), ^∗)≤Θ^∗-Θ^(k)_F/√(L_s T).

Furthermore, with the error bound for Θ^∗ - Θ^(k)_F shown in (<ref>), i.e. with probability at least (1 - ϵ) (1 - 15/a^2)^3s( 1 - 2/d - 14/m - 208/dm)^T, we have 

    sin∠ (^(k), ^∗)≤ O ( √(σ^2d/L_sm)).



The proof of Theorem <ref> is provided in Appendix <ref>.



  §.§.§ Remarks on the Main Theoretical Results

Theorem <ref> provides valuable insights into the error bound of Θ^(k) and Θ^∗. This result places significant reliance on the RIP condition, which in turn necessitates certain conditions to be met. Specifically, for the restricted isometry constant in Theorem <ref> to be small, it implies that the value of m should not be small (m≥Ω(s)). Additionally, the distance of the initial point should not be excessively large. In practice, the choice of the zero point is often suitable and aligns with the selection made in the subsequent numerical experiments. Importantly, this error bound diminishes as m increases.

In Theorem <ref>, we give a bound of the step size γ. Since we use a lot of reduction methods in the proof of this theorem, the available step size range can be larger in practice.

In practice, we often focus more on the distance metric sin∠(^(k),^∗). Theorem <ref> demonstrates that this distance can be bounded by the distance presented in Theorem <ref>, and it also tends to decrease as m increases. Notably, this distance should also decrease with an increase in T, as observed in the numerical experiments in next section. While it has been empirically observed, further theoretical validation remains a topic for future research.

Finally, we comment on how the theoretical results scale with sample size per task m and compare the scaling with that of several other algorithms. Recall that we require m = Ω(s log s), which results in a smaller m requirement compared to other methods. For instance, Burer-Monteiro factorization in <cit.> necessitates m≥Ω(s^4log(T)), and alternating minimization in <cit.> requires m≥Ω(s^2log(T)). The Method of Moments in <cit.> indeed shares a similar requirement with our approach, where it requires m = Ω (s log s). The error bound of the obtained Θ for this method, as specified in Theorem <ref>, is O(σ√(sσ^2sd+T/m)+s√(d/m)). This error bound may not be as favorable, as evident from the numerical results presented in the subsequent section. Another method, nuclear norm minimization in <cit.>, has been shown to require m≥Ω(1), which is favorable in theory but may not perform as well in practical scenarios, especially when s is small. The error bound of the obtained Θ for this method, as specified in Theorem <ref>, is O ( σ√(s ( d^2/m^2 + T/m))+√(s d/mmax{d, T}/m)). This bound is less favorable when the noise level σ is small. In summary, our method  holds a distinct advantage when dealing with scenarios where the sample complexity per task m is small.



§ NUMERICAL RESULTS





 §.§ Simulated Data



In this section, we conduct extensive numerical experiments to evaluate our proposed method, benchmarked against several competing methods, including both the classical model-based methods such as Method of Moments() and the modern model-agnostic methods such as . Specifically, we consider a setup where d = 100 and s = 5, following the experimental setup outlined in <cit.>, <cit.> and <cit.>. We use the following performance metric to compare different methods: (1) the normalized squared Frobenius norm distance

    𝐃𝐢𝐬𝐭_1 (Θ^†, Θ^∗) Θ^† - Θ^∗_F^2/T

and (2) the sine angle distance

    𝐃𝐢𝐬𝐭_2 (^†, ^∗) sin∠(^†, ^∗)

where Θ and , respectively, denote arbitrary values that Θ and  can take and Θ^∗ and ^∗ denote the ground truths, respectively. Using repeating draws from the true DGP, we investigate on average how these performance metrics are affected by key parameters including the task-specific sample size m, the total number of tasks T, and the noise level σ. Additionally, we examine their dynamics changing with iteration steps and system running time in a single experiment. More details on the experimental setups can be found in the Appendix.

We compare our method  with several alternative approaches:

    
  *  and : The alternating minimization and alternating minimization gradient descent algorithms studied in <cit.>.
    
  * : The ANIL algorithm outlined in <cit.>, which utilizes all available data in each iteration.
    
  * : The Burer-Monteiro factorization method as presented in <cit.>.
    
  * : The Method of Moments algorithm described in <cit.>.
    
  * : Algorithm 2 presented in <cit.>.
    
  * : The nuclear norm minimization method proposed in <cit.>.
    
  * : Our  method.






Figures <ref> and <ref> showcase how 𝐃𝐢𝐬𝐭_1 and 𝐃𝐢𝐬𝐭_2 change by varying the number of tasks T. When m = 25, it is notable that , , and  achieve lower errors compared to other methods. Notably, , designed to model-agnostically identify the meta-representation, consistently underperformed by other methods under the squared distances induced by the matrix Frobenius norm. When m = s = 5, the scenario corresponding to the extreme data-scarcity setting, methods like , , and  fail to learn either the regression coefficients (the left panel) or the task-invariant representation (the right panel). In this scenario, our proposed approach  exhibits superior performance, over other methods that still exhibit shrinking estimation error as T grows, including , , , and . 



Figure <ref> illustrates how the two metrics change as m varies. Notably, our proposed method  consistently maintains strong performance across different m's. Particularly, for small m values (m < 20), , , and  exhibit diminished performance. However, as m increases, the performance of  and  becomes better, and is eventually comparable to .



Figure <ref> illustrates how the two metrics change as σ changes. With m = 25, , , and  perform better than other methods in general. Additionally, these three methods, along with , exhibit the ability to reduce errors as σ diminishes. More comprehensive numerical results regarding situations where m is very small or σ=0 can be found in the Appendix <ref>. Notably,  fails to achieve zero error in the noiseless setting, since the implemented algorithm provided in <cit.> did not solve the optimization problem well[As explained in the Introduction, BFGS method is directly used to solve the optimization problem, which can lead to sub-optimal numerical performance. Theoretically, supported by the results of <cit.>,  should achieve zero error in the noiseless setting if the numerical optimization step achieves negligible iteration errors.].





Figure <ref> along with Table <ref> present the empirical minimum amount of data for each method to attain a sine angle distance ≤ 0.1 over one simulation. Notably,  and  requires an extensive amount of data and are thus excluded from this representation. Strikingly,  exhibits the most favorable performance in terms of data requirement, particularly when confronted with small sample size m per task.



Figure <ref> shows how 𝐃𝐢𝐬𝐭_1 and 𝐃𝐢𝐬𝐭_2 change over iterations and time for the four iterative methods: , ,  and . Evidently,  demonstrates significantly better computational efficiency compared to the other methods. Although  requires fewer iterations, its efficiency is hampered by the larger computation time per iteration. Notably, there are instances, as demonstrated in Appendix <ref>, where  still requires more iterations than  when dealing with small values of m.

In general,  exhibits the highest computational efficiency and consistently achieves superior statistical accuracy across various scenarios. While  and  are shown to be effective when data is not excessively scarce, they are ill-suited when the sample size of each task is small. Furthermore,  demands substantial computational resources per iteration. Its variant, , offers improved speed but lags behind other methods in terms of performance based on the two evaluation metrics. , , , and  can effectively handle scenarios with highly restricted data. Although  performs comparably to our method (), it falls short when confronted with extremely low noise intensity. On the other hand,  only provides a subspace estimation, yielding less satisfactory results. Despite the conciseness of  and , their demanding data requirements and suboptimal performance undermine their utility.



 §.§ Real Data Analysis



In this section, we further examine the performance of  and other competing methods in a real dataset. In particular, we use the air quality data from weather stations across different regions of China in 2023, which is originally from the China national urban air quality real-time publishing platform of the China Environmental Monitoring Station[The website of China national urban air quality real-time publishing platform is https://air.cnemc.cn:18007/https://air.cnemc.cn:18007/ and the dataset can be downloaded from https://quotsoft.net/air/#archivehttps://quotsoft.net/air/#archive.]. Model (<ref>) is adopted to study the relationship between the average PM2.5 in a day and the hourly CO value, hourly NO2 value, date, and geographical location. Specifically, the regression model for each weather station is treated as a single task, where the average value of 24-hour PM2.5 is the response variable (i.e., y), and the hourly CO value (24-dimensional), hourly NO2 value (24-dimensional), date (1-dimensional), and (geographical) location coordinates (2-dimensional) are combined into covariates or features (i.e., ), together with the intercept (all 1 vector). After removing corrupted data, outliers, and abnormal tasks, we are left with T = 1,210 tasks, each with a maximum sample size of 293 and a minimum sample size of 87 (i.e. m_t∈ [87, 293] for t = 1, ⋯, T). In order to better align with the modeling assumptions and the algorithmic implementation in this paper, we preprocess the data in the following steps:

    
  * Normalize the location coordinates of different tasks within the range [-1,1].
    
  * Sequentially assign date values x from 1,2,⋯,365 for each task. Since PM2.5 is generally low in summer and high in winter, we transform the date values to x=|x-183|, and then standardize them (x=x-x/σ_x, where x and σ_x are the mean and variance computed from the data).
    
  * Take the logarithm of the CO and NO2 values and standardize each dimension in every task.
    
  * Take the logarithm of PM2.5 values.

We assume that the coefficients of all tasks lie in an r-dimensional space. 80% of the tasks are used to train the model to get the low-dimensional subspace (), referred to as meta tasks. The remaining tasks are used to verify the validity of the resulting subspace, referred to as test tasks. For each meta task and test task, we divide the data into training points (randomly selecting m points) and test points (the rest). Training points in meta tasks are trained with different methods to obtain the subspace  and cross-task regression coefficients . Test points in meta tasks are used to evaluate the results of . Given , training points (_j,_j) in a test task are used to obtain the intrinsic low-dimensional coefficients by computing

    _j=((_j)^⊤(_j))^-1((_j)^⊤_j)

and

    _j = _j.

Then test points in test tasks are used to verify the performance of the obtained coefficients, which also reflects the performance of . The training and testing process is divided into four stages: use training points in meta tasks to train the model with different methods(meta-train); use test points in meta tasks to test the model(meta-test); and use training points in test tasks to solve the regression problem with (test-train); use test points in test tasks to test the performance(test-test).

We focus on the results of the stage meta-test, test-train and test-test. For each stage, we compute the mean relative error (MRE) of the predicted PM2.5 values for each task, and use the mean of MRE (M-MRE) over all corresponding tasks as the performance metric. Additionally, we use two other methods as baselines in test tasks:

    
  * :  is randomly generated when solving regression problems in test tasks;
    
  * : Combine the least square method and pseudo inverse, i.e. _j=(_j^⊤_j)^†_j^⊤_j.


We set r=6 and the results of different methods[ failed in the experiments of real data, so it is not shown in the results. Other approaches mentioned in Section <ref> are all compared here.] for various m are shown in Tables <ref>, <ref>, <ref>. The best-performing results are highlighted in red bold font and the second-best results are highlighted in blue bold font.







We also investigate the performance of the methods for different r values. The results are shown in Tables <ref>, <ref>, and <ref> for the meta-test, test-train, and test-test stages, respectively. In these tables, m is set to 30, and the best results are highlighted in red bold font, while the second-best results are in blue bold font.







The results in the tables show that our method  consistently ranks among the top two methods in different situations and is often the best one. Additionally, the methods based on Model (<ref>) significantly outperform the baseline methods, indicating that Model (<ref>) could be appropriate for investigating the given dataset.



§ CONCLUDING REMARKS
 


In this study, we introduced the concept of utilizing matrix rank minimization techniques to address challenges posed by multi-task linear regression problems with limited data availability. Our proposed approach, termed  (), demonstrates efficacy in resolving such problems. Moreover, we conducted a theoretical analysis of the error bounds associated with the estimations produced by this algorithm. Through experimental validation,  showcased superior accuracy and efficiency compared to existing methodologies, particularly in scenarios with exceptionally scarce data. This underscores the potential for improved solutions in cases involving very limited samples.

Nevertheless, it is important to acknowledge that, based on our experimental findings, the error bounds we established may not be the most tight. Notably, the error associated with the sine angle distance displays a gradual reduction with increasing values of T. This intriguing observation warrants further exploration and validation in future research endeavors.

In this paper, we assume that the rank s of the common subspace is known. The immediate next step is to show that our procedure adapts to unknown s. Another interesting direction is to analyze related algorithms, including , under the so-called “proportional asymptotic limit” <cit.> (s ≍ d ≍ n and s ≤ d), an arguably more realistic high-dimensional regime than the d ≫ n or d ≪ n, d →∞ regimes, using techniques from statistical physics <cit.>. 
It is also of theoretical interest to characterize the fundamental limit of estimating the common subspace when the task diversity assumption may be nearly violated, using the local minimax rate framework recently put forth in the theoretical computer science and statistics literature <cit.>. In most modern applications, the HPS linear model is at best an approximation of the actual data generating mechanism. Finally, developing methods for its non- or semi-parametric extension, such as HPS single index models or model-agnostic settings <cit.>, is another important research direction that deserves more attention.



§ ACKNOWLEDGMENTS

This work is supported by NSFC Grants  No.12090024 (CZ, XZ, LL) and No.12471274 (LL), Shanghai Municipal of Science and Technology Grants 21JC1402900 (LL) and 2021SHZDZX0102 (XZ, LL).

alpha





§ PROOFS
 


 §.§ Proof of Theorem <ref>
 

We analyze the concentration result for each term in Equation (<ref>). Using Lemma <ref>, we set _j to be the matrix 1/m(_t, j_t, j^⊤^⊤-_r). Consequently, =1/m_t^⊤_t^⊤-_r. 

We denote _l as the l-th row vector of  and _l_2=1. Since _t,j∼(0_d,_d), it follows that [(_l^⊤_t,j)^2]=1. According to [|z|^k]≤(2σ^2)^k/2kΓ(k/2) for any z∼(0,σ^2), we obtain [(_l^⊤_t,j)^2]≤15. By applying Chebyshev's Inequality, we can infer that ℙ(1-a ≤ (_l^⊤_t,j)^2≤ 1+a)≥1-15/a^2 for any a>0. Consequently, we have ^⊤_t,j_2^2=∑_l=1^r(_l^⊤_t,j)^2≤ (a+1)r with probability at least (1-15/a^2)^r for any a>0. 

Subsequently, _j can be bounded as follows:

    _j = 1/m(_t,j_t,j^⊤^⊤-)≤1/m_t,j_t,j^⊤^⊤ + 1/m = 1/m_t,j_2^2 + 1/m≤(a+1)r+1/m.

This holds with high probability. Since _j is Hermitian, we can compute the matrix variance statistic ν() under the assumptions stated in Lemma <ref>:

    ν() = ∑_j=1^m𝔼_j^2

The variance of each term is calculated as:

    𝔼_j^2     = 1/m^2( [_t,j(_t,j^⊤^⊤_t,j)_t,j^⊤^⊤]-) 
        = 1/m^2( [_t,j_2^2_t,j_t,j^⊤^⊤]-) 
       ≼1/m^2{(a+1)r [_t,j_t,j^⊤^⊤]-}
       ≼(a+1)r-1/m^2

Thus, 𝔼_j^2≤(a+1)r-1/m^2. As a result, we find that:

    ν() ≤∑_j=1^m𝔼_j^2 = (a+1)r-1/m

Applying Lemma <ref> with the aforementioned conditions and letting u=δ_r^(t)(assuming u≤1), we have

    1/m_t^⊤_t^⊤-≤δ_r^(t)≤√(8(a+1)r-4/3mlog2r/ϵ).

This holds with probability at least (1-ϵ)(1-15/a^2)^r for any a>0. If it holds for all i, we can set δ_r()=√(2(a+4)r+8/3mlog2r/ϵ) and sum up Equation (<ref>) over i. This leads to the conclusion that Equation (<ref>) holds with the probability at least (1-ϵ)^T(1-15/a^2)^rT, which in turn implies that Equation (<ref>) holds with this probability.




 §.§ Proof of Theorem <ref>
 
To start the proof, we first recall some relevant definitions and results from <cit.> and <cit.>.


Given a set of r rank-one matrices = {ψ_1,⋯,ψ_r}, there exists a set of s orthonormal matrices Γ = {γ_1, ⋯, γ_s} in the Frobenius sense[That is, ⟨γ_i, γ_j ⟩_ = 0 if i  j and γ_i_ = 1 for all i, such that (Γ) =  (). Here the inner product between two matrices is simply Frobenius inner product.]. We call Γ an orthonormal basis for the space (). We use P_ΓΘ to denote the projection of Θ onto the space (Γ). By definition, P_ΓΘ≡ P_Θ and
rank (P_ΓΘ) ≤ r, ∀Θ.



Assume that the rank-r matrix Θ_r has the singular value decomposition Θ_r=∑_i=1^r σ_iu_iv_i^⊤. Γ:=u_1v_1^⊤,u_2v_2^⊤,⋯,u_rv_r^⊤ is
called an SVD basis for the matrix Θ_r. Note that the elements in Γ are orthonormal rank-one matrices.



Suppose that the linear operator :ℝ^T × d→ℝ^N × 1 satisfies the RIP with constant δ_r(). Let  be an arbitrary orthonormal subset of ℝ^T × d such that
rank(P_Θ) ≤ r, ∀Θ∈ℝ^d × T. Then, for all b∈ℝ^N × 1 and Θ∈ℝ^T × d the following properties hold:

    P_^∗b≤√(1+δ_r())b_2, 
        (1-δ_r())P_Θ_≤P_^∗ P_Θ_≤ (1+δ_r())P_Θ_




Suppose that the linear operator :ℝ^T × d→ℝ^N × 1 satisfies the RIP with constant δ_r(). Let ,^' be an arbitrary orthonormal subset of ℝ^T × d such that
rank(P_∪^'Θ) ≤ r, ∀Θ∈ℝ^T × d. Then the following property hold:

    P_^∗(I-P_)Θ_≤δ_r()(I-P_)Θ_, ∀Θ∈ (^')




Suppose  is the best rank-s approximation to the matrix , and Γ is an SVD basis of . Then for any rank-s matrix _s and SVD basis Γ_s of _s, we have

    ‖ P_B-P_B‖_F ≤‖ P_B_s-P_B‖_F,

where B is any orthonormal set of matrices satisfying (Γ∪Γ_s)⊆ B.



    Since  is the best rank-s approximation to the matrix  and the rank of _s is s, ‖-‖_F ≤‖_s-‖_F. Hence,
    
    ‖ P_B(-) ‖_F^2 + ‖ (I-P_B)(-) ‖_F^2 ≤‖ P_B(_r-) ‖_F^2 + ‖ (I-P_B)(_r-) ‖_F^2.

    Since (I-P_B)=0 and (I-P_B)_r=0, this reduces to the result.


Provided here is a pivotal lemma that sets apart the proof in this paper from the approach presented in <cit.>.


Given the linear operator =𝒜/√(m) and the sub-Gaussian noise vector , the following inequality holds with probability at least (1-2/d-14/m-208/dm)^T:

    ^∗_≤^∗_≤√(2dTσ^2).




As ^∗_^2=∑_t=1^T_t^⊤_i_2^2/m, we can analyze _t^⊤_t_2^2/m for each task t individually. The expectation of it is [_t^⊤_t_2^2/m]= [1/m_t^⊤_t_t^⊤_t]=1/m [_t^⊤ [_t_t^⊤]_t]=dσ^2. According to [|z|^k]≤(2σ^2)^k/2kΓ(k/2) for any z∼(0,σ^2), the corresponding variance is 

    [_t^⊤_i_2^2/m]     =  [(∑_j=1^m∑_k=1^m_t,j^⊤_t,kε_t,jε_t,k)^2/m^2]-( [_t^⊤_t_2^2/m])^2 
       ≤ [15m(d^2+14d)+(m^2-m)(d^2+2d)]σ^4/m^2-d^2σ^4 
        = (2dm^2+14d^2m+208dm)σ^4/m^2.

By Chebyshev’s inequality, we can establish 

    ℙ(dσ^2-t≤_t^⊤_t_2^2/m≤ dσ^2+u)≥1-(2dm^2+14d^2m+208dm)σ^4/u^2m^2

for any u > 0. Here, We take u = dσ^2, resulting in _t^⊤_t_2^2/m≤ 2dσ^2 with a probability at least 1-2/d-14/m-208/dm. By summing up these inequalities over t, the lemma is proved.


Then we give the proof of Theorem <ref>.


Assuming that the Restricted Isometry Property (RIP) holds, as established in Theorem <ref>, let Γ^∗ and Γ^k be the SVD basis of Θ^∗ and Θ^(k), respectively. We denote an orthonormal basis of the subspace (Γ^∗∪Γ^k) as _k 
Starting with the relation Θ^∗-Θ^(k+1)_ = P__k+1Θ^∗-P__k+1Θ^(k+1)_ and Lemma <ref> (note that one can simply conduct a verbatim application of the lemma once replacing _r by ^*.), we can further expand this as follows:

    Θ^∗-Θ^(k+1)_
     
        =   P__k+1Θ^∗-P__k+1Θ^(k+1)_
    
        =   P__k+1Θ^∗-P__k+1Θ̂^(k+1)+P__k+1Θ̂^(k+1)-P__k+1Θ^(k+1)_
    ≤   P__k+1Θ^∗-P__k+1Θ̂^(k+1)_ + P__k+1Θ̂^(k+1)-P__k+1Θ^(k+1)_
    ≤    2P__k+1Θ^∗-P__k+1Θ̂^(k+1)_

Here, we define Z^k=Θ^∗-Θ^(k), which results in Θ̂^(k+1)=Θ^(k)+γ^∗( Z^k+/√(m)). Continuing the derivation, we obtain:

    P__k+1Θ^∗-P__k+1Θ̂^(k+1)_
    
        =    P__k+1Θ^∗-P__k+1Θ^(k)-γ P__k+1^∗ Z^k-γ P__k+1^∗/√(m)_
    
        =    P__k+1Z^k-γ P__k+1^∗ P__k+1Z^k-γ P__k+1^∗(I-P__k+1)Z^k-γ P__k+1^∗/√(m)_
    ≤   (I-γ P__k+1^∗ P__k+1)P__k+1Z^k_ + γP__k+1^∗(I-P__k+1)Z^k_ + γ/√(m)P__k+1^∗_.

By using γ≤1 and utilizing Proposition <ref> and <ref>, the following inequalities are established with a probability at least (1-ϵ)(1-15/a^2)^3s(1-2/d-14/m-208/dm)^T for any a>0:





    P__k+1Θ^∗-P__k+1Θ̂^(k+1)_
    ≤    (1-γ+γδ_2s()) P__k+1Z^k_ + γδ_3s()(I-P__k+1)Z^k_ + γ/√(m)^∗_
    ≤    (1-γ+γδ_3s()) P__k+1Z^k_ + (1-γ+γδ_3s())(I-P__k+1)Z^k_ + γ/√(m)^∗_
    ≤   √(2)(1-γ+γδ_3s())Z^k_ + √(2dTσ^2/m).

This implies Z^k+1_F≤ 2√(2)(1-γ+γδ_3s())Z^k_ + √(8dTσ^2/m). Consequently, with 1-1/2√(2)/1-δ_3s<γ<1, i.e. 0<2√(2)(1-γ+γδ_3s())<1, we have:

    Θ^∗-Θ^(k)_
    ≤   [2√(2)(1-γ+γδ_3s())]^kZ^0_ + {[2√(2)(1-γ+γδ_3s())]+[2√(2)(1-γ+γδ_3s())]^2+⋯
        +[2√(2)(1-γ+γδ_3s())]^k}√(8dTσ^2/m)
    ≤   [2√(2)(1-γ+γδ_3s())]^kZ^0_ + 2√(2)/2√(2)(1-γ+γδ_3s())√(dTσ^2/m)

In particular, when γ=1, with the result in Theorem <ref>,

    Θ^∗-Θ^(k)_
    ≤    (2√(2)δ_3s())^kZ^0_ + 2√(2)/2√(2)δ_3s()√(dTσ^2/m)
    ≤   (192(a+1)s-32/3mlog6s/ϵ)^k/2Z^0_ + 2√(2)σ/1-√(192(a+1)s-32/3mlog6s/ϵ)√(dT/m).


If k is large enough, i.e. with k→∞, we have

    Θ^∗-Θ^(k)_≤ O(√(σ^2dT/m)).





 §.§ Proof of Theorem <ref>
 
We will draw on a similar result, namely Lemma 15 in <cit.>, and adopt their methodologies to prove Theorem <ref>.

We begin the proof by defining f ()  - ^∗^∗_^2 for short. Since ^∗ minimizes f (), it follows from the first-order condition f/=0 that ^∗=^∗^∗^⊤. Therefore for any , we have ^∗^∗^⊤-^∗^∗_^2≤-^∗^∗_^2. Moreover,

    ^∗^∗^⊤-^∗^∗_^2    = ^∗^∗(^⊤-)_^2 
       = ^∗^∗_^⊤__^2 
       = ^∗^∗_^⊤_^2 
       = tr(^∗^⊤^∗^∗_^⊤_^∗^⊤) 
       ≥λ̅_s(^∗^⊤^∗)^∗_^⊤_^2 
       ≥λ̅_s(^∗^∗^⊤)^∗_^⊤_2^2.

In the above inequalities, λ̅_s (A) represents the s-th largest eigenvalue of a given matrix A. It is worth noting that the second-to-last inequality follows from the property that for any positive semi-definite matrices M and N, (MN)≥σ_min(M) (N). Let Θ^(k) and ^(k) be the output obtained by  after k iterations. As λ̅_s(^∗^∗^⊤)=λ_sT, it follows that:

    sin∠(^(k),^∗) = ^∗_k^⊤_2 ≤Θ^∗-Θ^(k)_/√(λ_s T)≤Θ^∗-Θ^(k)_/√(L_s T).

Combining this result with (<ref>), we immediately have (<ref>).







§ MORE EXPERIMENTAL RESULTS AND DETAILS OF SIMULATED DATA EXPERIMENTS
 
In this part, we show more numerical results and give the details of all the experiments in Section <ref>.


 §.§ More Experimental Results
 
We previously demonstrated the trend of metrics' evolution with varying T for cases where m=25 and m=5, depicted in Figures <ref> and <ref>. Now, we extend our analysis to the scenario where m=10, as illustrated in Figure <ref>. Similar to our previous observations, our proposed method  continues to outperform others under these conditions. Moreover, , , and  exhibit strong performance when T is relatively large. This underscores the importance of having both sufficiently large values of m and T to harness the effectiveness of these three methods.



In Figure <ref>, we present results identical to those shown in Figure <ref>, with the sole exception of setting σ=0. Notably, , , , and  achieve remarkably low error rates, aligning with our expectations. In contrast,  underperforms in this scenario, while , , and  exhibit even poorer performance. These findings are in line with the outcomes depicted in Figure <ref>.



We also investigated a more challenging scenario with s=25. Figure <ref> and Figure <ref> illustrate the evolution of the metrics with respect to T for cases when m=40 and m=25. The comparison of results closely resembles that of the s=5 scenario with m=10 and m=5. These findings suggest that the methods can effectively address such challenges even as s increases.





Figure <ref> displays the same experiment as shown in Figure <ref>, except for T=1600. The results in these two scenarios are notably similar.



The above experiments investigate the influence of both m and T on the results. These two parameters determine the amount of data utilized to address the problems. All of these experiments consistently demonstrate that our method, , requires less data to solve the same problem and consistently achieves superior results under the same conditions.

We have presented results pertaining to the parameter σ in Figure <ref>. Figure <ref> illustrates the scenario when m=10 and T=1600, while Figure <ref> covers the case when m=5 and T=6400. Notably, , , and  are omitted in Figure <ref> because they cannot effectively function when m=5. In Figure <ref>, the results closely resemble those in Figure <ref>, with the exception that  exhibits less effectiveness. This observation aligns with the findings in Figures <ref> and <ref>. The outcomes presented in Figure <ref> underscore that our method, , stands out as the sole approach capable of performing effectively in scenarios with extremely limited data represented by small values of m.





Figure <ref> serves as an illustrative example of how the metrics evolve with iterations and time. To provide further insights, we present three additional examples under the same conditions of m=25 and T=400 in Figure <ref>. Furthermore, in Figure <ref>, we replicate the experiments with m=10 and T=3200. Notably,  exhibits suboptimal performance in the second example, signifying its instability when dealing with small values of m. Collectively, these experiments underscore that our method, , not only delivers effective results but also achieves efficiency by necessitating fewer iterations and less time to yield favorable outcomes.







 §.§ Experimental Details
 
In all the experiments, we generate the true task-invariant matrix ^∗ by first QR factorizing of a d × d matrix with elements sampled i.i.d. from the standard normal 𝒩 (0, 1), and then retrieving the first s columns. The elements of each task-varying coefficient vector _i^∗ are also sampled i.i.d. from the standard normal 𝒩 (0,1). For all methods except , initialization is based on random draws from the standard normal as well. In contrast,  initializes _i^(0) as the zero vector for all i ∈ [T]. For the algorithm , a hyperparameter λ called the “regularization coefficient” is determined as λ = σ/T√(T + d^2 / m/mT), following the recommendation in <cit.>. When implementing the  and  algorithms, as in <cit.> and <cit.>, we use the L-BFGS algorithm <cit.> within the scipy package. Moreover, we employ the autograd package for gradient computations in these two methods.

The choice of step sizes and maximum iterations for our proposed , as well as , , and , are determined empirically. These parameters are set differently for various experiments. In Figure <ref>, the step sizes for  and  are 0.25 and 1.0, respectively, while the inner and outer loop step sizes for  are both 0.5. , ,  and  are run for 40, 200, 20, 200 iterations, respectively.
In Figure <ref>, Figure <ref>, Figure <ref> and Figure <ref>, the step sizes for  and  are 0.05 and 0.2, respectively, while the inner and outer loop step sizes for  are both 0.1. , ,  and  are run for 200, 1000, 100, 800 iterations, respectively.
In Figure <ref>, the step sizes for  and  are 0.5 and 1.0, respectively, while the inner and outer loop step sizes for  are both 0.5. , ,  and  are run for 300, 200, 100, 400 iterations, respectively.
In Figure <ref> and <ref>, the step sizes for  and  are 0.5 and 1.0, respectively, while  the inner and outer loop step sizes for  are both 0.5. , ,  and  are run for 20, 50, 20, 50 iterations, respectively.
In Figure <ref> and Figure <ref>, the step sizes for  and  are 0.1 and 0.4, respectively, while  the inner and outer loop step sizes for  are both 0.2. , ,  and  are run for 100, 500, 50, 400 iterations, respectively.
In Figure <ref>, the step sizes for  and  are 0.25 and 1.0, respectively, while  the inner and outer loop step sizes for  are both 0.5. , ,  and  are run for 1000, 200, 100, 400 iterations, respectively. In addition, the regularization coefficient for  is set to 1 × 10^-5 since the best value is zero when it follows the formulation in <cit.>.
In Figure <ref>, the step sizes for  and  are 0.3 and 0.4, respectively, while  the inner and outer loop step sizes for  are both 0.2. , ,  and  are run for 15000, 500, 400, 500 iterations, respectively.
In Figure <ref>, the step sizes for  and  are 0.2 and 0.4, respectively, while  the inner and outer loop step sizes for  are both 0.3. , ,  and  are run for 40, 100, 40, 100 iterations, respectively.

To enhance the reliability of our results, all the results in Figure <ref> to Figure <ref> and Figure <ref> to Figure <ref> are averaged over 5 trials. The results presented in Figure <ref> and Table <ref> are derived from experimenting with different T values for different m, approximated to the nearest hundred or ten, and represent the average error over 5 trials for each scenario.





