additions

# Multi-Task Reinforcement Learning in Continuous Control with Successor
Feature-Based Concurrent Composition

Yu Tang Liu1,2 and Aamir Ahmad2,1 1Max Planck Institute for Intelligent
Systems, 72076 Tübingen, Germany. 2Institute of Flight Mechanics and Controls,
University of Stuttgart, 70569 Stuttgart, Germany.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Deep reinforcement learning (DRL) frameworks are increasingly used to solve
high-dimensional continuous-control tasks in robotics. However, due to the
lack of sample efficiency, applying DRL for online learning is still
practically infeasible in the robotics domain. One reason is that DRL agents
do not leverage the solution of previous tasks for new tasks. Recent work on
multi-tasking DRL agents based on successor features has proven to be quite
promising in increasing sample efficiency. In this work, we present a new
approach that unifies two prior multi-task RL frameworks, SF-GPI and value
composition, for the continuous control domain. We exploit compositional
properties of successor features to compose a policy distribution from a set
of primitives without training any new policy. Lastly, to demonstrate the
multi-tasking mechanism, we present a new benchmark for multi-task continuous
control environment based on Raisim. This also facilitates large-scale
parallelization to accelerate the experiments. Our experimental results in the
Pointmass environment show that our multi-task agent has single task
performance on par with soft actor critic (SAC) and the agent can successfully
transfer to new unseen tasks where SAC fails. We provide our code as open-
source at https://github.com/robot-perception-group/concurrent_composition for
the benefit of the community.

##  I Introduction

One approach to achieve optimal continuous control in robotics is through
reinforcement learning-based methods (RL) [1]. However, RL agents are usually
optimized for one specific task and require more training if the reward
definition changes. On the other hand, multi-task RL frameworks aim to design
an RL agent that allows recycling old policies for other tasks to achieve
higher sample efficiency[2].

One way to achieve multi-task RL is through transfer learning by manipulating
the policies trained on the old tasks. In this work, we focus on the transfer
learning methods that define tasks using linearly decomposable reward
functions, commonly used to define many robot control tasks [3].

Figure 1: In concurrent composition, policy extraction is too expensive to
perform in every time step. Instead, we propose composing the primitives
directly in the run time.

Currently, there are two major transfer learning frameworks following the
aforementioned assumption. The first framework, using successor features and
generalized policy improvement (SF-GPI, [4, 5]), derives the new policy by
choosing the sub-policies, or primitives, with the largest associated action-
value. The other framework directly composes the action-value of different
tasks, which we refer to as value composition (VC [6, 7, 8, 9]). This
framework first defines the composite value functions of new tasks by linearly
combining the constituent action-value functions and then extracts the policy.

Different from the Option framework[10], which is an another approach to
multi-task RL that combines the primitives in temporal order, in our novel
approach we focus on concurrent policy composition by composing the primitives
at every time step. The samples collected by the composite policy expose each
primitive to a higher quality data and thus enable the knowledge transfer and
increase the sample efficiency [4]. Policy extraction is prohibitive since it
requires updating policies multiple times. We consider real-time performance
as critical. Therefore, we must adapt VC-based methods for online composition.
We do so by multiplicative policy composition (MCP [11]). We show that this
can be achieved via the successor features framework, where SF-GPI becomes one
special case.

In summary our contributions include the following.

  * •

Derivation of the relation between value composition and policy composition.

  * •

A novel method unifying the GPI and VC under the successor feature-based
online concurrent composition framework.

  * •

Introduction of a new benchmark environment based on Raisim for multi-task
reinforcement learning.

##  II Related work

The compositional agents in our framework build upon the two transfer learning
frameworks: VC and SF-GPI. VC-based methods mainly came from the maximum
entropy framework[12]. By encouraging policy randomness, the resulting policy
is usually more robust against the noise [13] and exhibits better exploratory
behavior[14]. The value composition methods often exploit the property that
the optimal value function of conjunction tasks can be well approximated by
the average of the constituent value functions [6, 7]. Interestingly, boolean
operators can further expand the valid compositions in the discrete setting
and achieve a super-exponential amount of task transfer[8, 9]. However, none
of these approaches consider concurrent composition scenarios.

On the other hand, SF-GPI[4] applies successor features [15] to estimate a set
of Q-functions. Maximizing these Q-functions results in the policy
improvement. The results in [4] indicate that training several primitives
simultaneously can significantly improve the sample efficiency more than
training individually. Our work can be seen as an extension to SF-GPI by
bringing valid operations from VC to successor feature-enabled concurrent
learning.

##  III Background

###  III-A Multi-Task Reinforcement Learning

RL improves a control policy through interaction with the environment, and it
formulates this procedure by Markov decision processes (MDPs),
$\mathcal{M}\equiv<\mathcal{S},\mathcal{A},p,r,\gamma>$. $\mathcal{S}$ and
$\mathcal{A}$ define the state space and action space respectively. The
dynamics $p(s_{t+1}|s_{t},a_{t})$ defines the transition probability from one
state to another. Reward function $r(s_{t},a_{t},s_{t+1})$ governs the desired
agent behavior. The discount $\gamma\in[0,1)$ determines the valid time
horizon of the task. The goal is to find the optimal control policy
$\pi(a|s):S\rightarrow A$, such that the expected discounted return,
$J(\pi)=\mathbb{E}_{\pi,\mathcal{M}}[\sum_{\tau=t}^{\infty}\gamma^{\tau-t}r_{\tau}]$,
can be maximized.

We consider transfer within a set of MDPs $\mathcal{M}_{w}\in\mathcal{T}$
where $\mathcal{M}_{w}\equiv<\mathcal{S},\mathcal{A},p,r_{w},\gamma>$ , and
each has its unique reward function,
$r_{w}(s_{t},a_{t},s_{t+1})=\phi(s_{t},a_{t},s_{t+1})^{\top}\cdot w$, which is
a linear combination of a set of commonly shared $d$ features,
$\phi(s_{t},a_{t},s_{t+1})\in\mathbb{R}^{d}$ and the task weight
$w\in\mathbb{R}^{d}$. The goal is to find a set of primitives that can solve
the subset of tasks $\mathcal{T}^{\prime}\in\mathcal{T}$ such that we can
recover the optimal policies from the unseen tasks by the found primitives.

###  III-B Maximum Entropy RL

The maximum RL framework includes a bonus policy entropy term in the expected
reward function. The goal is to find an optimal control policy that not only
maximizes the reward but also encourages the policy randomness,

$J(\pi)=\mathbb{E}_{\pi,\mathcal{M}_{w}}\left[\sum_{\tau=t}^{\infty}\gamma^{\tau-t}(r_{w,\tau}+\alpha\mathcal{H}(\pi(\cdot|s_{\tau})))\right],$
(1)

where $\alpha$ is a temperature scaler controlling the entropy of the policy
and $\mathcal{H}(\cdot)$ is an entropy estimator. Prior works have found that
the bonus term can encourage exploration [16], increase the training policy’s
robustness against noise [13], and preserve the policy multimodality [14].

One tractable way to achieve maximum entropy RL is through temporal difference
learning by boostrapping the value estimation. Here, we define the action-
value function, or Q-function, of policy $\pi$ on task $w$ as:

$Q^{\pi}_{w}(s_{t},a_{t})\equiv
r_{w,t}+\mathbb{E}_{\pi}\left[\sum_{\tau=t+1}^{\infty}\gamma^{\tau-t}(r_{w,\tau}+\alpha\mathcal{H}(\pi(\cdot|s_{\tau})))\right].$
(2)

Note when $\alpha=0$, we recover the standard RL Q-function. The Q-function
can be perceived as a performance metric of its policy and the estimation
quality can be improved recursively through the soft policy iteration.

###  III-C Soft Policy Iteration

The soft policy iteration can be defined in two steps: (1) soft policy
evaluation (2) soft policy improvement [14]. The soft policy evaluation step
allows to improve the quality of Q-value estimation and it can be formulated
as follows,

$\displaystyle Q^{\pi}_{w}(s_{t},a_{t})$
$\displaystyle=r_{w}(s_{t},a_{t})+\gamma\mathbb{E}_{s_{t+1}\sim
p}\left[V(s_{t+1})\right]$ (3) $\displaystyle V^{\pi}_{w}(s_{t})$
$\displaystyle=\mathbb{E}_{a_{t}\sim\pi}\left[Q^{\pi}_{w}(s_{t},a_{t})-\log\pi(a_{t}|s_{t})\right].$
(4)

The policy distribution can be updated according to the soft policy
improvement step:

$\pi_{new}=\operatorname*{\arg\\!\min}_{\pi^{\prime}}D_{KL}\left(\pi^{\prime}(\cdot|s_{t})||\frac{\exp(Q^{\pi_{old}}_{w}(s_{t},\cdot))}{Z^{\pi_{old}}(s_{t})}\right),$
(5)

where $D_{KL}$ denotes the Kullback-Leiber divergence and the normalizing
factor, $Z^{\pi_{old}}(s_{t})$, can be neglected because it does not
contribute to the gradient computation. Intuitively, the goal is to match the
policy distribution to the Q-value estimation of different actions; therefore,
the higher the Q-value, the higher the probability of selecting the action,
i.e.

$\pi(\cdot|s)\propto e^{Q^{\pi}(s,\cdot)}.$ (6)

Alternating between the policy evaluation step and policy improvement step
will eventually converge to the optimal Q-function $Q^{*}_{w}$ on task $w$
with its corresponding optimal policy $\pi^{*}$, where $Q^{*}_{w}\equiv
Q^{\pi^{*}}_{w}$ [14].

###  III-D Successor Feature (SF)

Given that
$r_{w}(s_{t},a_{t},s_{t+1})=\bm{\phi}(s_{t},a_{t},s_{t+1})^{\top}\cdot w$, we
can rewrite the Q-function in following form,

$\displaystyle Q^{\pi}_{w}($ $\displaystyle
s_{t},a_{t})=\mathbb{E}_{\pi}\left[\sum_{\tau=t}^{\infty}\gamma^{\tau-t}r_{w,\tau}\right]$
(7)
$\displaystyle=\mathbb{E}_{\pi}\left[\sum_{\tau=t}^{\infty}\gamma^{\tau-t}\bm{\phi}^{\top}_{t}\right]\cdot
w=\bm{\psi}^{\pi}(s_{t},a_{t})^{\top}\cdot w$ (8)

where $\bm{\psi}^{\pi}\in\mathbb{R}^{d}$ is the successor feature of policy
$\pi$[15]. To approximate it, we can utilize function approximators and apply
soft value iteration by replacing the reward with features in analogous to
(3). We introduce the successor features to the maximum entropy framework by
including the entropy term. As a result, we obtain the policy iteration for
the successor features,

$\displaystyle\bm{\psi}^{\pi}(s_{t},a_{t})$
$\displaystyle=\bm{\phi}_{t}+\mathbb{E}_{\pi}\left[\sum_{\tau=t}^{\infty}\gamma^{\tau-t}(\bm{\phi}_{\tau}+\alpha\mathcal{H}(\pi(\cdot|s_{\tau}))\right]$
(9) $\displaystyle=\bm{\phi}_{t}+\gamma\mathbb{E}_{s_{t+1}\sim
p(\cdot|s_{t},a_{t})}\left[\bm{\Upsilon}(s_{t+1})\right],$ (10)
$\displaystyle\bm{\Upsilon}^{\pi}(s_{t})$
$\displaystyle=\mathbb{E}_{a_{t}\sim\pi}\left[\bm{\psi}^{\pi}(s_{t},a_{t})-\log\pi(a_{t}|s_{t})\right].$
(11)

The policy improvement step is similar to (5) by replacing the Q-function with
the approximation from the successor features and the desired primitive task
weight, which is defined in Sec. IV-B.

###  III-E Generalized Policy Improvement Composition (GPI)

Given policies $\Pi=\\{\pi_{1},\pi_{2},...,\pi_{n}\\}$ and their corresponding
Q-functions in a new task $w$,
$Q^{\Pi}_{w}=\\{Q^{\pi_{1}}_{w},Q^{\pi_{2}}_{w},...,Q^{\pi_{n}}_{w}\\}$ . If a
new policy $\pi^{GPI}$ is constructed in the following way,

$\pi^{GPI}(s)\in\operatorname*{\arg\\!\max}_{a\in\mathcal{A}}\max_{\pi\in\Pi}Q^{\pi}_{w}(s,a),$
(12)

then the GPI theorem [4] states that the new policy $\pi^{GPI}$ is no worse
than all other policies since its Q-function is always greater than its
members, i.e. $Q^{\pi^{GPI}}_{w}(s,a)\geq\sup_{\pi\in\Pi}Q^{\pi}_{w}(s,a).$

This composition is not the optimal composition and the composite value is an
under-estimate of the true value
$Q^{GPI}_{w}=\max_{\pi\in\Pi}Q^{\pi}_{w}(s,a)\leq Q^{*}_{w}(s,a)$. And
subsequently, the composite policy performance might not be optimal
$Q^{\pi^{GPI}}_{w}\leq Q^{\pi^{*}}_{w}$ [7].

Successor features provide one quick and convenient way to evaluate
performance of a policy $\pi$ in any task $w$. Therefore, SF-GPI algorithm can
evaluate all policies’ performance on the new task $w$ instantaneously by the
SFs, called generalized policy evaluation, or GPE. Then GPI can be applied to
select the action that leads to the highest value.

Note the $\operatorname*{\arg\\!\max}_{a\in\mathcal{A}}$ expression is not
tractable in continuous action space and we instead sample actions from the
primitives, i.e. $\operatorname*{\arg\\!\max}_{a\sim\Pi}$. We describe our
implementation of the continuous SF-GPI in Sec. IV-C3.

###  III-F Value Composition

Combining the Q-functions allows to obtain Q-function of the new task without
solving them at all. Prior works have defined a set of valid operations in
multi-goal discrete control [8, 9]. In continuous control, however, we focus
on the linear task combinations which do not assume binary task weights. That
is, if a new task is a linear combination of old tasks, $r_{new}=a\cdot
r_{a}+b\cdot r_{b}$, then the Q-value of the new task can be approximated by
$\hat{Q}^{\sum}_{new}\approx a\cdot Q^{*}_{a}+b\cdot Q^{*}_{b}$. The new
policy can be extracted from the composite Q-function. However, since policy
extraction would require several policy improvement steps, it is intractable
during interaction with the environment.

Note that the composition is an overestimate of the true optimal value
function, i.e. $\hat{Q}^{\sum}_{new}\geq Q^{*}_{new}$. Since it assume all
constituent tasks can be achieved at the same time across the whole state
space [7, 6]. When this is not the case, the composite policy’s is not
optimal, i.e. $Q^{\pi^{\sum}}_{new}\leq Q^{\pi^{*}}_{new}$.

###  III-G Multiplicative Compositional Policy (MCP)

MCP [11] has shown to compose the policy primitive in the immitation learning
scenario. Here, we introduce it to concurrent reinforcement learning. Given a
set of $d$ Gaussian primitives $\Pi\in\mathcal{N}^{dim(\mathcal{A})\times d}$
and a gating vector $g\in[0,1]^{d}$ where $||g||_{1}=1$, the MCP function
$MCP(\Pi,g):\mathcal{N}^{dim(\mathcal{A})\times
d}\times\mathbb{R}^{d}\rightarrow\mathcal{N}^{dim(\mathcal{A})\times 1}$ is
defined as follows,

$\pi(a|s)=MCP(\Pi,g)\equiv\frac{1}{Z(s)}\prod_{i=0}^{d}\pi_{i}(a|s)^{g_{i}}.$
(13)

where $Z(s)$ is a partition function to normalize the distribution. If all the
primitives are Gaussian, then the composite is also a Gaussian. The author did
not explain why this composition is valid, but we show that this is a natural
outcome from value composition in the policy space (Sec. IV-C). For
convenience, we refer to the MCP-based compositional method as compositional
policy improvement, or CPI, in contrast to GPI.

##  IV Methodology

Our novel approach is described in the following sub-sections.

###  IV-A Successor Feature-based Composition

In contract to the GPI composition (Sec. IV-C3), we introduce two more SF-
based compositions.

####  IV-A1 Successor Feature based Value Composition (SFV)

First we introduce the linear value composition (Sec. III-F) to the successor
feature framework. Given the task $r=\bm{\phi}^{\top}\cdot w$, we can treat
features $\bm{\phi}$ as sub-tasks. From this perspective, $\bm{\phi}$ becomes
the constituent sub-tasks and $r$ is the composite task. Thus, the successor
features then represent the policy’s performance on each sub-task. Therefore,
(8) can be seen as a special case for value composition of the same policy.

Since we do not know the task policy $\pi$, we can only approximate it by the
existing primitives $\bm{\Pi}=[\pi_{0},...,\pi_{n}]$. Value composition
theorems state that the average of the Q-function is a sufficiently good
approximation of the true Q-function [6, 8, 9].

$\displaystyle Q^{\pi}_{w}(s,a)\approx
Q^{{\sum}}_{w}(s,a)=\frac{1}{n}\sum_{i=0}^{n}Q_{w}^{\pi_{i}}(s,a)=\frac{\bm{w}^{\top}}{n}\sum_{i=0}^{n}\bm{\psi}^{\pi_{i}}(s,a)$
(14)

From this expression, one can see that SFV composition can only achieve an
average performance in all sub-tasks. Consider optimal SF, since SFV assumes
all tasks can be achieved at the same time, this is also an overestimate of
the true value function, i.e. $\hat{Q}^{SFV}_{w}\equiv\hat{Q}^{\sum}_{w}\geq
Q^{*}_{w}$. We can define the greedy policy extracted by the following
composition.

$\pi^{SFV}(s)\in\operatorname*{\arg\\!\max}_{a\in\mathcal{A}}\frac{\bm{w}^{\top}}{n}\sum_{i=0}^{n}\bm{\psi}^{\pi_{i}}(s,a).$
(15)

Note when the sub-tasks are not simultaneously achievable, this policy can
only achieve a sub-optimal performance $Q^{\pi^{SFV}}_{w}(s,a)\leq
Q^{\pi^{*}}_{w}(s,a)$.

####  IV-A2 Maximum Successor Feature Composition (MSF)

If we can treat successor features as the sub-tasks, then it naturally raises
the question: can we combine primitives in a certain way such that we achieve
the current best achievable performance in all sub-tasks and subsequently any
composite task? That is, given SFs of their corresponding $n$ primitives,
$\Psi=\\{\bm{\psi}^{\pi_{0}},\bm{\psi}^{\pi_{1}},...,\bm{\psi}^{\pi_{n}}\\}$,
where
$\bm{\psi}^{\pi_{i}}=\\{\psi_{0}^{\pi_{i}},\psi_{1}^{\pi_{i}},...\psi_{d}^{\pi_{i}}\\}$.
Then the desired composition operator $\mathcal{U}$ can be defined as
$\mathcal{U}\psi^{\pi}\equiv\\{\max_{\pi\in\Pi}{\psi_{0}^{\pi}},...,\max_{\pi\in\Pi}{\psi_{d}^{\pi}}\\}$.
With this, we define a composition rule,

$\pi^{MSF}(s)\in\operatorname*{\arg\\!\max}_{a\in\mathcal{A}}w^{\top}\mathcal{U}\psi^{\pi}(s,a).$
(16)

Considering optimal SF, then this composition will have the largest over-
estimation than any other compositions, i.e.

$\hat{Q}^{MSF}_{w}\geq\hat{Q}^{SFV}_{w}\geq Q^{*}_{w}\geq Q^{GPI}_{w}.$ (17)

###  IV-B Construct Policy Primitives

We construct policy primitives based on the sub-tasks and each primitive is
only responsible for its own sub-task. Therefore, we have the same number of
sub-tasks and primitives (i.e. $d=n$). All primitives are trained by soft
policy iterations and follow its convergence guarantee, i.e.

$\displaystyle\pi_{i}(\cdot|s)\propto
e^{Q^{\pi_{i}}_{\mathbf{e}_{i}}(s,\cdot)}=e^{\mathbf{e}_{i}^{\top}\cdot\psi^{\pi_{i}}(s,\cdot)}.$
(18)

The base tasks $\mathbf{e}_{i}\in\mathbb{R}^{d}$ can be constructed
arbitrarily, as long as the base task and the transfer task share the same
sub-tasks. For example, the base task defined as the zero vector will fail
since this task has no overlapping with the transfer task. For simplicity, we
define it as a zero vector with the i-th component equal to 1. Our goal is to
compose these primitives such that the true task policy can be recovered.

###  IV-C Approximate True Policy by Composing Primitives

This is the core part of our method. We provide analytical view on how the
task policy can be reconstructed by a set of task primitives given their
relation in value space introduced in Sec. IV-A.

####  IV-C1 MSF Composition

Following MSF in value space (16) we derive the correspondent policy
composition,

$\displaystyle\pi(\cdot|s)$ $\displaystyle\propto
e^{Q^{\pi}_{w}(s,\cdot)}=e^{w^{\top}\cdot\bm{\psi}^{\pi}(s,\cdot)}\approx
e^{w^{\top}\cdot\mathcal{U}\psi^{\pi}(s,a)}$ (19)
$\displaystyle=e^{w^{\top}\cdot[\psi^{\pi_{0}}_{0},\psi^{\pi_{1}}_{1},...\psi^{\pi_{d}}_{d}](s,\cdot)}\propto\prod_{i=0}^{d}\pi_{i}(\cdot|s)^{w_{i}}.$
(20)

The approximation in (19) assumes the maximum SF of each sub-task is most
likely to be the SF of the primitive trained from the respective sub-task.
Then the MSF composition can be achieved by MCP, i.e.
$\hat{\pi}(\cdot|s)=MCP(\Pi,g^{MSF})$, where $g^{MSF}=\frac{w}{|w|}$ is the
normalized gating vector.

Algorithm 1 MSF Composition

1:function MSF(s, w):

2: $\hat{\pi}(\cdot|s)=MCP(\Pi(\cdot|s),\frac{w}{|w|})$ $\triangleright$ CPI

3: $a\sim\hat{\pi}(\cdot|s)$ $\triangleright$ sample from composite

4: return $a$

####  IV-C2 SFV Composition

Following SFV in value space (14) we derive the correspondent policy
composition,

$\displaystyle\pi(\cdot|s)\propto e^{Q^{\pi}_{w}(s,\cdot)}\approx
e^{Q^{{\sum}}_{w}(s,\cdot)}=\prod_{i=0}^{d}e^{\frac{\bm{w}^{\top}}{d}\bm{\psi}^{\pi_{i}}(s,\cdot)}$
(21)
$\displaystyle=\prod_{i=0}^{d}e^{\frac{1}{d}w_{i}\psi_{i}^{\pi_{i}}(s,\cdot)}\prod_{i=0}^{d}\prod_{j\neq
i}^{d}e^{\frac{1}{d}w_{j}\psi_{j}^{\pi_{i}}(s,\cdot)}\propto\prod_{i=0}^{d}\pi_{i}(\cdot|s)^{\frac{w_{i}}{d}}\cdot\epsilon$
(22)

where $\epsilon$ is the artefactual noise. Because of the noise, the SFV
composition is not a valid composition in policy space with the set of
primitives introduced in Sec. IV-B.

####  IV-C3 GPI Composition

Though this is already a SF-based composition, we reformulate it for
completeness.

$\displaystyle\pi(\cdot|s)\propto e^{Q^{\pi}_{w}(s,\cdot)}$
$\displaystyle\approx e^{\max_{\pi_{i}\in\Pi}Q^{\pi_{i}}_{w}(s,\cdot)}$ (23)
$\displaystyle\propto\operatorname*{\arg\\!\max}_{\pi_{i}\in\Pi}Q^{\pi_{i}}_{w}(s,\cdot)$
(24)

The approximation can be achieved by successor features, i.e.
$\hat{\pi}(\cdot|s)=\operatorname*{\arg\\!\max}_{\pi_{i}\in\Pi}Q^{\pi_{i}}_{w}(s,\cdot)=\operatorname*{\arg\\!\max}_{\pi_{i}\in\Pi}w^{\top}\psi^{\pi_{i}}(s,\cdot)$.

Algorithm 2 GPI Composition

1:function SF-GPI(s, w):

2: $\bm{a}\sim\Pi(\cdot|s)$ $\triangleright$ sample from primitives

3: $\hat{\bm{Q}}^{\pi}_{w}=\bm{\psi}^{\pi}(s,\bm{a})\cdot w$ $\triangleright$
GPE

4: $a=\operatorname*{\arg\\!\max}_{a\sim\Pi}\hat{\bm{Q}}^{\pi}_{w}$
$\triangleright$ GPI

5: return $a$

At this point, we have constructed an analytic method to find the policy
composition corresponding to the composition in value space. Note that if
primitives are constructed with different base tasks $\mathbf{e}_{i}$, the
composition in policy space will change accordingly. In addition, the VC-based
method will always end up with MCP composition in policy space, and therefore,
they fall into the CPI category. Next. we heuristically extend the idea to
direct action composition, which allows composition happens at action
component level.

###  IV-D Multiplicative Compositional Action (MCA)

Rewriting the Gaussian primitives $\pi_{i}\in\Pi$ as the following,

$\displaystyle\pi_{i}(\cdot|s)=[\mathcal{N}^{\pi_{i}}_{0},...,\mathcal{N}^{\pi_{i}}_{\mathcal{A}}|s]$
(25)

where $\mathcal{N}^{\pi_{i}}_{j}=\mathcal{N}^{\pi_{i}}(\mu_{j},\sigma_{j})$ is
the shorthand for the j-th action component of $\pi_{i}$ . Then the MSF
composition (20) has the following form,

$\displaystyle\hat{\pi}(\cdot|s)=MCP(\Pi,g)=[\prod_{i=0}^{d}{(\mathcal{N}^{\pi_{i}}_{0})^{g_{i}}},...,\prod_{i=0}^{d}{(\mathcal{N}^{\pi_{i}}_{\mathcal{A}})^{g_{i}}}].$
(26)

However, since all action components are coupled by the scale $g_{i}$, we
cannot manipulate them differerntly. Alternatively, given
$\kappa\in[0,1]^{dim(\mathcal{A})\times n}$ where $||\kappa||_{1}=1$, we
define the direct action composition function
$MCA(\Pi,\kappa):\mathcal{N}^{dim(\mathcal{A})\times
n}\times\mathbb{R}^{dim(\mathcal{A})\times
n}\rightarrow\mathcal{N}^{dim(\mathcal{A})\times 1}$:

$\displaystyle\hat{\pi}(\cdot|s)=MCA(\Pi,\kappa)\equiv[\prod_{i=0}^{d}{(\mathcal{N}^{\pi_{i}}_{0})^{\kappa^{\pi_{i}}_{0}}},...,\prod_{i=0}^{d}{(\mathcal{N}^{\pi_{i}}_{\mathcal{A}})^{\kappa^{\pi_{i}}_{\mathcal{A}}}}].$
(27)

We define $\kappa^{\pi_{i}}_{j}$ in such a way that action components with
lower successor feature value will be filtered out and higher value to be
emphasized, i.e. $\kappa^{\pi_{i}}\propto\mathcal{P}_{\mathcal{A}\times
d}\cdot\bm{\psi}^{\pi_{i}}$, where $\mathcal{P}$ is an impact matrix defining
the mapping from successor feature to the action components.

###  IV-E Impact Matrix

We would like to know how each action component is related to each feature
throughout the whole trajectory, that is, the successor feature. Then given
some successor features $\bm{\psi}=\\{\psi_{0},\psi_{1},...\psi_{d}\\}$, we
can define the impact matrix as absolute value of a Jacobian of $\psi$ w.r.t
action, i.e.

$\displaystyle\mathcal{P}\equiv|\nabla_{a}{\psi(s,a)}|=\begin{bmatrix}|\frac{\partial\psi_{0}}{\partial
a_{0}}|&...&|\frac{\partial\psi_{d}}{\partial a_{0}}|\\\ \vdots&\ddots&\\\
|\frac{\partial\psi_{0}}{\partial
a_{\mathcal{A}}}|&...&|\frac{\partial\psi_{d}}{\partial
a_{\mathcal{A}}}|\end{bmatrix}_{\mathcal{A}\times d}$ (28)

The major concern of this expression is the amount of computation it requires
and accuracy of the derivative. Luckily, in practice, this estimation can be
computed efficiently by vectorized computation via GPUs. Yet the accuracy
depends mostly on the function smoothness. Fortunately, the value estimation
usually becomes smoother during training and increasing the robustness of the
function approximator can create similar effect as well [17]. But still, the
estimation can be very noisy in the early phase of training. We can attenuate
the noise by averaging over all policies
$\bm{\psi}\approx\frac{1}{n}\sum_{n}{\bm{\psi}^{\pi_{n}}}$. Therefore, the
more primitives the less noisy it will become. To further alleviate the noise,
we include the temporal dependency to the estimation by averaging with the
impact matrix from the previous time step, i.e.
$\hat{\mathcal{P}}\approx(\mathcal{P}_{t-1}+\mathcal{P}_{t})/2$.

###  IV-F Direct Action Composition (DAC)

With the impact matrix, we define $\kappa\equiv\mathcal{P}_{\mathcal{A}\times
d}\cdot\Psi_{d\times n}\odot w_{d\times 1}$, where

$\displaystyle\Psi=\\{\bm{\psi}^{\pi_{0}},\bm{\psi}^{\pi_{1}},...,\bm{\psi}^{\pi_{n}}\\}=\begin{bmatrix}\psi_{0}^{\pi_{0}}&...&\psi_{0}^{\pi_{n}}\\\
\vdots&\ddots&\\\
\psi_{d}^{\pi_{0}}&...&\psi_{d}^{\pi_{n}}\end{bmatrix}_{d\times n},$ (29)

$\odot$ denotes the Hadamard product, and $n=d$. Intuitively, task weights $w$
first filter the task-relevant successor features and then map their impact to
each action components by the impact matrix. In practice, we replace $\Psi$
with advantage $\Gamma$ for numerical stability, i.e.
$\Gamma=\Psi-\frac{1}{n}\sum_{n}{\bm{\psi}^{\pi_{n}}}$. To make sure
$||\kappa||_{1}=1$, $softmax$ is applied. Empirically, we found the
performance much better when we filter out the advantage below zero, i.e.
$\Gamma_{clip}=max(0,\Gamma)$. Finally we have

$\displaystyle\hat{\kappa}=softmax(\mathcal{P}\cdot\Gamma_{clip}\odot w),$
(30)

then we can achieve DAC by $MCA(\Pi,\hat{\kappa})$ as defined in (27).

Algorithm 3 DAC Compostion

1:function DAC(s, w):

2: $\bm{a}\sim\Pi(\cdot|s)$ $\triangleright$ sample from primitives

3:
$\hat{\kappa}=softmax(\hat{\mathcal{P}}(s,\bm{a})\cdot\Gamma_{clip}(s,\bm{a})\odot
w)$

4: $\hat{\pi}(\cdot|s)=MCA(\Pi(\cdot|s),\hat{\kappa})$ $\triangleright$ CPI

5: $a\sim\hat{\pi}(\cdot|s)$ $\triangleright$ sample from composite

6: return $a$

Lastly, we introduce another variant, DAC-GPI, by choosing the action
component with the largest $\hat{\kappa}$, i.e.

$\displaystyle\hat{\pi}(\cdot|s)\equiv[\mathcal{N}^{\operatorname*{\arg\\!\max}_{\pi\in\Pi}{\hat{\kappa}^{\pi}_{0}}}_{0},...,\mathcal{N}^{\operatorname*{\arg\\!\max}_{\pi\in\Pi}{\hat{\kappa}^{\pi}_{\mathcal{A}}}}_{\mathcal{A}}].$
(31)

Algorithm 4 DAC-GPI Compostion

1:function DAC-GPI(s, w):

2: $\bm{a}\sim\Pi(\cdot|s)$ $\triangleright$ sample from primitives

3:
$\hat{\kappa}=softmax(\hat{\mathcal{P}}(s,\bm{a})\cdot\Gamma_{clip}(s,\bm{a})\odot
w)$

4:
$a=[\operatorname*{\arg\\!\max}_{a_{0}}\hat{\kappa}_{0}^{\pi},...,\operatorname*{\arg\\!\max}_{a_{\mathcal{A}}}\hat{\kappa}_{\mathcal{A}}^{\pi}]$
$\triangleright$ GPI

5: return $a$

Notice that when all the action components come from the same policy, DAC-GPI
reduces to SF-GPI. Moreover, when each primitive shares the same exponent
across all action components, DAC reduces to MSF. The major difference between
DAC and DAC-GPI is that the former compose action components from all
primitives according to the value while the latter only selects the best one.
Therefore, we expect DAC to scale better with increasing amount of features.

In summary, we have constructed a method to bridge the composition in value
space to policy space. This provides an effective analytic method to generate
different kinds of concurrent composition algorithms. We then show that the
composition can not only happen in policy space but also in action space with
the help of the impact matrix. Different from standard RL agents, the
compositional agents usually involve critics or SFs to the action selection
process and filter out the lower quality actions before delivery.

###  IV-G Training the Network

The training data is gathered by direct interaction with the environment by
the concurrent composition. Therefore, we can treat the concurrent
compositional agent just as any regular RL agent and train in an end-to-end
fashion without any pre-training.

In total, we have $n$ pairs of SFs and primitives with the $i$-th successor
feature network $\psi^{\pi_{i}}$, target network $\bar{\psi}^{\pi_{i}}$, and
the primitive network $\pi_{i}$ paramterized by $\theta_{\psi^{\pi_{i}}}$,
$\theta_{\bar{\psi}^{\pi_{i}}}$, and $\theta_{\pi_{i}}$ respectively. The
training scheme follows the soft policy iteration, which includes the
evaluation step (10,11) and improvement step (5). The evaluation step now
evaluates the successor feature instead of Q-function while the policy
improvement step stays the same by following the base task Q-function
$Q^{\pi_{i}}_{\mathbf{e}_{i}}$. Therefore, the loss for the i-th successor
feature can be computed as follows,

$\displaystyle J_{\bm{\psi}^{\pi_{i}}}(\theta_{\bm{\psi}^{\pi_{i}}})$
$\displaystyle=\mathbb{E}_{(s_{t},a_{t})\sim\mathcal{D}}\left[(\bm{\psi}^{\pi_{i}}(s_{t},a_{t})-\hat{\bm{\psi}}^{\pi_{i}}(s_{t},a_{t}))^{2}\right],$
(32)

where $\mathcal{D}$ is the replay buffer and we have the TD-target
$\hat{\bm{\psi}}^{\pi_{i}}$ approximated by the target network
$\bar{\psi}^{\pi_{i}}$,

$\displaystyle\hat{\bm{\psi}}^{\pi_{i}}(s_{t},a_{t})$
$\displaystyle=\bm{\phi}_{t}+\gamma\mathbb{E}_{s_{t+1}\sim
p(\cdot|s_{t},a_{t})}\left[\bm{\Upsilon}^{\pi_{i}}(s_{t+1})\right],$ (33)
$\displaystyle\bm{\Upsilon}^{\pi_{i}}(s_{t})$
$\displaystyle=\mathbb{E}_{a_{t}\sim\pi_{i}}\left[\bm{\bar{\psi}}^{\pi_{i}}(s_{t},a_{t})-\log\pi_{i}(a_{t}|s_{t})\right].$
(34)

Then the primitive can be improved by reducing the KL-divergence in (5), which
can be approximated by minimizing the loss,

$\displaystyle
J_{\pi_{i}}(\theta_{\pi_{i}})=\mathbb{E}_{s_{t}\sim\mathcal{D},a_{t}\sim\pi_{i}}\left[\log\pi_{i}(a_{t}|s_{t})-Q^{\pi_{i}}_{\mathbf{e_{i}}}(s_{t},a_{t})\right],$
(35)

with the primitive action value defined by (18).

To reduce the effect of off-policyness, we apply retrace [18], a variant of
importance sampling with less variance, in SFs estimation. By assuming that
each primitive has the best performance in its own sub-task, the correction
can be performed conveniently by comparing the action distribution with other
primitives instead of the behavior policy, which has complicated compositional
distribution.

We summarize our novel approach by Algorithm. 5. The method alternates between
the interaction phase to collect samples by the compositional policy and
learning phase that updates the function approximators’ parameters by
stochasitc gradient descend with samples collected previously. In the
algorithm, looping over all primitives seems intimidating, but in practice we
can vectorize the computation and update all primitives simultaneously in one
pass. Therefore, the total amount of computation only scales linearly with the
amount of primitives and the computation time will be the same with a GPU.

Algorithm 5 Composition Agent

1:Choose $method\in$ MSF, SF-GPI, DAC, DAC-GPI

2:Initialize network parameters
$\bm{\theta_{\psi^{\pi}}},\bm{\theta_{\bar{\psi}^{\pi}}},\bm{\theta_{\pi}}$

3:for N steps do

4: if exploration then $a\leftarrow$ Uniform($\mathcal{A}$)

5: else$~{}a\leftarrow method(s)$

6:
$\mathcal{D}\leftarrow\mathcal{D}\cup(s_{t},a_{t},r(s_{t},a_{t}),\phi(s_{t},a_{t}),s_{t+1})$

7: for each gradient step do

8: for $i\leftarrow 1,2,...,n$ do

9:
$\theta_{\psi^{\pi_{i}}}\leftarrow\theta_{\psi^{\pi_{i}}}-\lambda_{\psi}\nabla_{\theta_{\psi^{\pi_{i}}}}J_{\psi^{\pi_{i}}}(\theta_{\psi^{\pi_{i}}})$

10:
$\theta_{\pi_{i}}\leftarrow\theta_{\pi_{i}}-\lambda_{\pi}\nabla_{\theta_{\pi_{i}}}J_{\pi_{i}}(\theta_{\pi_{i}})$

11:
$\theta_{\bar{\psi}^{\pi_{i}}}\leftarrow\tau\theta_{\bar{\psi}^{\pi_{i}}}+(1-\tau)\theta_{\psi^{\pi_{i}}}$

lastly, for clarity, we note that in the proposed approach there exist 3 folds
of transfer mechanism and each complement with one another to accelerate the
learning. First, the SFs allows transfer to any arbitrary task by specifying
different task weights $w$. The task transfer enable the policy evaluation in
a new task without any training. This kind of transfer is associated with
composition in the value space. Second, since the composite policy is
guaranteed to have better performance than the primitives, the trajectory
collected by this composite is expected to explore higher value region in the
composite task. In this case, the knowledge is transferred from the expert
policy to the primitives, similar to teacher-student [19] or learn from
demonstration [20]. Third, if the primitives and SFs share representations,
such as in Fig. 2, the knowledge is transferred among themselves, and each SF
treat others as auxiliary tasks to improve the representations [21, 22]. Or,
given a new SF, it is likely that the learning will be faster than training
from scratch due to the pre-trained feature extraction layers.

##  V Experiments and Results

The experiment aims to understand the transfer of the primitives. Then we test
if the multi-task RL agents generated by the proposed framework, i.e., SF-GPI,
MSF, DAC, DAC-GPI, can be task-agnostic and transfer to any unseen task and
whether they can leverage knowledge transfer and speed up learning. We perform
all experiments on a single computer (AMD Ryzen Threadripper 3960X, 24x
3.8GHz, NVIDIA GeForce RTX2080 Ti, 11GB). Agents are implemented based on
Pytorch[23]. The hyper-parameter tuning tools are provided by Wandb[24].

(a) primitive network.

(b) successor feature network.

Figure 2: network architecture.

(a) The first row corresponds to the two action components from the primitive
$\pi_{x}=(a_{x},a_{y})$ and the second row for $\pi_{y}$ in the 2D plane. The
color bar represents the action magnitude between $[-1,1]$.

(b) The first row corresponds to the two components from kappa
$\kappa^{\pi_{x}}=(\kappa^{\pi_{x}}_{x},\kappa^{\pi_{x}}_{y})$ and the second
row for $\kappa^{\pi_{y}}$ in the 2D plane. The magnitude $[0,1]$ is
represented by the color bar.

Figure 3: Composition visualization of Pointmass2D-Simple

(a) The SF-GPI agent with different feature sets. The higher slope of the
learning curve implies a faster learning speed.

(b) When the evaluation task is different from the training one, the baseline
SAC[16] fails to generalize to the new task while multi-task agents are less
affected.

(c) The hyper-parameter tuning is conducted with the composite task
$w^{navigation}=(\mathbf{1},0.5,1)$ and $w^{hover}=(\mathbf{1},1,1)$ and
during evaluation, the success weight is increased to 100 for observing the
transfer ability. Here we show the best five sets of hyper-parameters.

(d) We compare composition agents with the baseline in the single task
$w=(\mathbf{1},0.5,1)$. Baselines SAC[16] is the SOTA maximum entropy RL agent
while SQL[14] applies an energy-based training method to achieve expressive
multi-modal policy.

Figure 4: Every curve is represented by the mean with two standard deviations
of 5 experiments run by the best agent from the hyper-parameter tuning. The
hyper-parameter search is conducted by bayesian optimization. We perform in
total a hundred searches for every agent. The ablation study indicates that
dropout layer can hurt the performance while the effect of other techniques
remain unclear, including layernorm, entropy tuning, importance sampling,
activation function, etc. Since DAC-GPI behaves fairly similarly to DAC and
MSF performs relatively poorly, we neglect them for some experiments.

###  V-A Network Architecture

We represent primitives and successor features by two neural networks as
displayed in Fig. 2. Since the layers or representations are shared, the
knowledge transfer can happen within the network. To reduce over-estimation in
the SFs and enable higher update-to-data ratio (UTD ratio) [25], we include
dropout and layernorm to the SF network. Since ReLU activation can result in
negative transfer [22], we incorporate SeLU[26] instead. To further mitigate
the positive bias [27], we apply the minimum double critic technique for the
SF estimation, i.e.,
$\hat{\psi}^{\pi}=min(\psi_{\theta_{1}}^{\pi},\psi_{\theta_{2}}^{\pi})$.

###  V-B Environment and Tasks

We build our training environment in the Raisim[28] simulator, which allows
massive vectorization for the environment to collect samples in parallel. Even
though Raisim has several built-in robots, to address our problem, we
introduce an abstract omnidirectional robot, Pointmass.

####  V-B1 Pointmass

We first introduce the state and action space of the environment,

  * •

state space $\mathcal{S}$:
$\Delta{\mathbb{X}},\Delta{\mathbb{V}},\Delta{\Theta},\Delta{\Omega}\in\mathbb{R}^{3}$,
are relative position, velocity, orientation, and angular velocity between the
robot and goal respectively.

  * •

action space $\mathcal{A}$: $(a_{x},a_{y},a_{z})$, are the applied force on
x-, y-, and z-axis of the Pointmass respectively.

Features, or sub-tasks, can be defined arbitrarily as long as they can be
derived from the state space. We pre-define three sets of features for our
experiments, as shown in Table. 1. The default will be Regular if not
explicitly mentioned.

Regular |  $(-|\Delta{\mathbb{X}}|,-||\Delta{\mathbb{V}}||_{2},Success)$
---|---
Simple |  $(-|\Delta{\mathbb{X}}|)$
Augment |  $(-|\Delta{\mathbb{X}}|,-||\Delta{\mathbb{X}}||_{2},-|\Delta{\mathbb{V}}|,-||\Delta{\mathbb{V}}||_{2},Success)$

TABLE 1: Feature sets for Pointmass environment. $|\cdot|$ is the absolute
value whereas $||\cdot||_{2}$ is the L2 norm. $Success$ is a binary feature
triggered when relative distance is less than 1 meter, i.e., $Success=1$ if
$||\mathbb{X}||_{2}<1$ else $0$.

The Regular feature set represents a standard way to define reward functions
in robotics, Simple is simplified for the demonstration purpose, and Augment
includes the redundant features, and each generates an extra redundant
primitive, e.g., both $-|\Delta{\mathbb{X}}|$ and
$-||\Delta{\mathbb{X}}||_{2}$ can achieve position control. The goal of each
primitive is to maximize its corresponding feature or sub-task. The effect of
features is displayed in Fig. 4(a). By adding more features, we obtain more
primitives with diverse behavior and increase the learning speed (the slope of
the curve).

Though the overall environment is defined in 3D space, it can be reduced to a
lower dimension by reducing the action space, e.g., Pointmass2D has the action
space $\mathcal{A}=(a_{x},a_{y})$, and in this case, the primitive associated
with altitude feature is redundant.

####  V-B2 Composite Task

We introduce a composite task, which includes two tasks during training,
navigation and hover task. The navigation task aims to bring the point mass
closer to the goal position, while the hover task requires the robot’s
velocity to approach zero when near the goal. The task switches from
navigation to hover when the distance to the goal is less than five meters. We
can formulate this accordingly,

$\begin{array}[]{l}w=\begin{cases}w^{navigation}&\text{if
$d(s_{pointmass},s_{target})\geq 5m$ }\\\
w^{hover}&\text{otherwise}\end{cases}\par\end{array}$ (36)

where $w^{navigation}=(\mathbf{1},0.5,1)$ and $w^{hover}=(\mathbf{1},20,1)$
and bold symbols denote the vector of the same length as the corresponding
feature. Because the velocity penalty is significant when the agent is within
the hovering range, the agent will stop immediately for both multi-tasking and
standard agents. Now, to test the agents’ transfer capability, in evaluation,
we change the composite task to emphasize the success reward, i.e.,
$w^{navigation}_{eval}=(\mathbf{1},0.5,10000)$ and
$w^{hover}_{eval}=(\mathbf{1},20,10000)$. We can demonstrate in Fig. 4(b) that
the multi-task agent is less affected when the evaluation task changes.

###  V-C Training Stability and Composition Noise

The composition noise stems from the action components redundant to the
transfer task and it could undermine the training stability. In
Pointmass2D-Simple, the primitives corresponding to the
features,$-|\Delta{X}|,-|\Delta{Y}|$ produce the distributions
$\pi_{x}(a|s)=(\mathcal{N}^{\pi_{x}}_{x},\epsilon)$ and
$\pi_{y}=(\epsilon,\mathcal{N}^{\pi_{y}}_{y})$, where $\epsilon$ is the
composition noise of unknown distribution (Fig. 3(a)). With MSF, the composite
contains the noise in both action components,
$\pi^{MSF}=(\mathcal{N}^{\pi_{x}}_{x}\epsilon,\mathcal{N}^{\pi_{y}}_{y}\epsilon)$,
whereas SF-GPI contains noise in one of the action components
$\pi^{GPI}=(\mathcal{N}^{\pi_{x}}_{x}\epsilon,\mathcal{N}^{\pi_{y}}_{y})$ or
$(\mathcal{N}^{\pi_{x}}_{x},\mathcal{N}^{\pi_{y}}_{y}\epsilon)$. Unlike MSF,
GPI can utilize critics to select the less noisy one.

Ultimately, only DAC and DAC-GPI can attenuate the noise in both action
components,
$\pi^{DAC}=((\mathcal{N}^{\pi_{x}}_{x})^{\kappa^{\pi_{x}}_{x}}(\epsilon)^{\kappa^{\pi_{y}}_{x}},(\mathcal{N}^{\pi_{y}}_{y})^{\kappa^{\pi_{x}}_{y}}(\epsilon)^{\kappa^{\pi_{y}}_{y}})$,
when the critics reduce $\kappa$ of the noisy term (Fig. 3(b) displays
$\kappa$). Fig. 4(c) reflects that the ability to filter out the composition
noise is critical to learning speed. DAC-GPI has the best convergence speed,
followed by DAC and then SF-GPI. We found the MSF composition unstable due to
the composition noise.

Lastly, compared to the maximum entropy RL baselines, i.e., SAC, SQL, in the
single task setting (Fig. 4(d)), the multi-task agents show a similar
convergence property as SAC.

##  VI Conclusions and Future Work

Our work unifies the SF-GPI and value composition to the continuous concurrent
composition framework and allows reconstructing task policy from a set of
primitives. The proposed method was extended to composition at the action
component level. We demonstrate in the Pointmass environment that our multi-
task agents can reconstruct the task policy from a set of primitives in real
time and transfer the skills to solve unseen tasks while the single-task
performance is competitive with SAC. This flexible framework incorporates well
with the reward-shaping techniques, such as entropy regularization,
curiosity[29], etc. In addition, the task-agnostic property should benefit the
autotelic framework [30] where agents can set goals and curriculum for
themselves [31].

However, the primary concern at this stage is whether the proposed approach
can scale to higher dimensional problems. Additionally, two important topics
are left as future works. First, look for the corresponding value composition
for DAC. A good starting point might be thinking of the MSF composition with
weights evaluated by GPE. Second, the optimality of each composition method.
One might start with bounding the loss incurred by the policy and value
composition.

## References

  * [1] J. Kober, J. A. Bagnell, and J. Peters, “Reinforcement learning in robotics: A survey,” _The International Journal of Robotics Research_ , vol. 32, no. 11, pp. 1238–1274, 2013.
  * [2] M. E. Taylor and P. Stone, “Transfer learning for reinforcement learning domains: A survey.” _Journal of Machine Learning Research_ , vol. 10, no. 7, 2009.
  * [3] T. Haarnoja, S. Ha, A. Zhou, J. Tan, G. Tucker, and S. Levine, “Learning to walk via deep reinforcement learning,” _Robotics: Science and Systems_ , 2019\.
  * [4] A. Barreto, S. Hou, D. Borsa, D. Silver, and D. Precup, “Fast reinforcement learning with generalized policy updates,” _Proceedings of the National Academy of Sciences_ , vol. 117, no. 48, pp. 30 079–30 087, 2020.
  * [5] A. Barreto, D. Borsa, S. Hou, G. Comanici, E. Aygün, P. Hamel, D. Toyama, S. Mourad, D. Silver, D. Precup _et al._ , “The option keyboard: Combining skills in reinforcement learning,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [6] T. Haarnoja, V. Pong, A. Zhou, M. Dalal, P. Abbeel, and S. Levine, “Composable deep reinforcement learning for robotic manipulation,” in _2018 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2018, pp. 6244–6251.
  * [7] J. Hunt, A. Barreto, T. Lillicrap, and N. Heess, “Composing entropic policies using divergence correction,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 2911–2920.
  * [8] B. Van Niekerk, S. James, A. Earle, and B. Rosman, “Composing value functions in reinforcement learning,” in _International conference on machine learning_. PMLR, 2019, pp. 6401–6409.
  * [9] G. Nangue Tasse, S. James, and B. Rosman, “A boolean task algebra for reinforcement learning,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 9497–9507, 2020.
  * [10] R. S. Sutton, D. Precup, and S. Singh, “Between mdps and semi-mdps: A framework for temporal abstraction in reinforcement learning,” _Artificial intelligence_ , vol. 112, no. 1-2, pp. 181–211, 1999.
  * [11] X. B. Peng, M. Chang, G. Zhang, P. Abbeel, and S. Levine, “Mcp: Learning composable hierarchical control with multiplicative compositional policies,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [12] B. D. Ziebart, _Modeling purposeful adaptive behavior with the principle of maximum causal entropy_. Carnegie Mellon University, 2010.
  * [13] R. Fox, A. Pakman, and N. Tishby, “Taming the noise in reinforcement learning via soft updates,” _Conference on Uncertainty in Artificial Intelligence. AUAI Press_ , 2015.
  * [14] T. Haarnoja, H. Tang, P. Abbeel, and S. Levine, “Reinforcement learning with deep energy-based policies,” in _International conference on machine learning_. PMLR, 2017, pp. 1352–1361.
  * [15] P. Dayan, “Improving generalization for temporal difference learning: The successor representation,” _Neural computation_ , vol. 5, no. 4, pp. 613–624, 1993.
  * [16] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _International conference on machine learning_. PMLR, 2018, pp. 1861–1870.
  * [17] P. Pauli, A. Koch, J. Berberich, P. Kohler, and F. Allgöwer, “Training robust neural networks using lipschitz bounds,” _IEEE Control Systems Letters_ , vol. 6, pp. 121–126, 2021.
  * [18] R. Munos, T. Stepleton, A. Harutyunyan, and M. Bellemare, “Safe and efficient off-policy reinforcement learning,” _Advances in neural information processing systems_ , vol. 29, 2016.
  * [19] S. Schmitt, J. J. Hudson, A. Zidek, S. Osindero, C. Doersch, W. M. Czarnecki, J. Z. Leibo, H. Kuttler, A. Zisserman, K. Simonyan _et al._ , “Kickstarting deep reinforcement learning,” _arXiv preprint arXiv:1803.03835_ , 2018.
  * [20] T. Hester, M. Vecerik, O. Pietquin, M. Lanctot, T. Schaul, B. Piot, D. Horgan, J. Quan, A. Sendonaris, I. Osband _et al._ , “Deep q-learning from demonstrations,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 32, no. 1, 2018.
  * [21] M. Jaderberg, V. Mnih, W. M. Czarnecki, T. Schaul, J. Z. Leibo, D. Silver, and K. Kavukcuoglu, “Reinforcement learning with unsupervised auxiliary tasks,” _arXiv preprint arXiv:1611.05397_ , 2016.
  * [22] H. Wang, E. Miahi, M. White, M. C. Machado, Z. Abbas, R. Kumaraswamy, V. Liu, and A. White, “Investigating the properties of neural network representations in reinforcement learning,” _arXiv preprint arXiv:2203.15955_ , 2022.
  * [23] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [24] L. Biewald, “Experiment tracking with weights and biases,” 2020, software available from wandb.com.
  * [25] T. Hiraoka, T. Imagawa, T. Hashimoto, T. Onishi, and Y. Tsuruoka, “Dropout q-functions for doubly efficient reinforcement learning,” _International Conference on Learning Representations (ICLR)_ , 2022.
  * [26] G. Klambauer, T. Unterthiner, A. Mayr, and S. Hochreiter, “Self-normalizing neural networks,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [27] S. Fujimoto, H. Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _International conference on machine learning_. PMLR, 2018, pp. 1587–1596.
  * [28] J. Hwangbo, J. Lee, and M. Hutter, “Per-contact iteration method for solving contact dynamics,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 2, pp. 895–902, 2018. [Online]. Available: www.raisim.com
  * [29] D. Pathak, P. Agrawal, A. A. Efros, and T. Darrell, “Curiosity-driven exploration by self-supervised prediction,” in _Proceedings of the 34th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, D. Precup and Y. W. Teh, Eds., vol. 70. PMLR, 06–11 Aug 2017, pp. 2778–2787.
  * [30] C. Colas, T. Karch, O. Sigaud, and P.-Y. Oudeyer, “Autotelic agents with intrinsically motivated goal-conditioned reinforcement learning: a short survey,” _Journal of Artificial Intelligence Research_ , vol. 74, pp. 1159–1199, 2022.
  * [31] S. Narvekar, B. Peng, M. Leonetti, J. Sinapov, M. E. Taylor, and P. Stone, “Curriculum learning for reinforcement learning domains: A framework and survey,” _The Journal of Machine Learning Research_ , vol. 21, no. 1, pp. 7382–7431, 2020.


# Fast randomized algorithms for computing the generalized tensor SVD based on
the tubal product

Salman Ahmadi-Asl  Ugochukwu Ugwu

###### Abstract

This work deals with developing two fast randomized algorithms for computing
the generalized tensor singular value decomposition (GTSVD) based on the tubal
product (t-product). The random projection method is utilized to compute the
important actions of the underlying data tensors and use them to get small
sketches of the original data tensors, which are easier to handle. Due to the
small size of the sketch tensors, deterministic approaches are applied to them
to compute their GTSVDs. Then, from the GTSVD of the small sketch tensors, the
GTSVD of the original large-scale data tensors is recovered. Some experiments
are conducted to show the effectiveness of the proposed approaches.

###### keywords:

Randomized algorithms, generalized tensor SVD, tubal product

###### MSC:

15A69 , 46N40 , 15A23

[label1]organization=Center for Artificial Intelligence Technology, Skolkovo
Institute of Science and Technology, Moscow, Russia<EMAIL_ADDRESS>[label2]organization=Department of Mathematics, Colorado State University,
Fort Collins, USA<EMAIL_ADDRESS>
##  1 Introduction

The Singular Value Decomposition (SVD) is one of the key matrix
factorizations, which has been widely used in many applications such as signal
processing and machine learning [1]. It can compute the best low-rank
approximation of a matrix in the least-squares sense for any invariant matrix
norm. This has been a perfect property convincing its wide applications. The
SVD is applied to a single matrix and can capture orthonormal bases for its
four fundamental subspaces. The idea of generalizing the SVD to a pair of
matrices was first proposed in [2, 3] and is called generalized SVD (GSVD)
with applications in solving inverse problems [4] and genetics [5, 6]. The
classical SVD or GSVD is prohibitive for computing low-rank approximations of
large-scale data matrices. The randomization approach is highly effective at
solving this issue when computing the SVD or GSVD of big matrices [7, 8, 9].
The principal idea of the randomization framework is first capturing the range
of a given data matrix by either multiplication with random matrices or
sampling some actual columns of the original data matrix. Then, using an
orthonormal basis for the mentioned range, a small sketch of the matrix is
computed, which can be processed and analyzed easily. Finally, the SVD/GSVD of
the original big data matrices is recovered from the SVD/GSVD of the small
sketches. This idea leads to a significant speed-up in our computations
although it also offers parallel implementation. It is known that randomized
algorithms can provide stable approximations and this issue together with the
other benefits mentioned above, have made them ubiquitous tools in numerical
linear algebra. Due to the promising results achieved by the randomized
algorithms for matrices, they have been extended to tensors and for different
types of tensor decompositions, see [10, 11, 12, 13].

In this work, we focus on the tensor SVD (t-SVD) proposed in [14, 15]. The
t-SVD has similar properties as the classical SVD. In particular, in contrast
to the Tucker decomposition [16, 17] or the canonical polyadic decomposition
[18, 19], its truncation version provides the best tubal rank approximation in
the least-squares sense. The GSVD was generalized to tensors based on the
t-SVD in [20, 21], which we call GTSVD. Motivated by the works [8, 9], we
develop two fast randomized algorithms for computing the GTSVD.

The key contributions of this work are as follows

  * 1.

Developing two fast randomized algorithms for the computation of the GTSVD
based on the t-product.

  * 2.

Providing computer simulations convincing the applicability of the proposed
randomized algorithm.

The rest of the paper is structured as follows. In Section 2, the
preliminaries are presented. Here, we also introduce the t-SVD model and the
necessary algorithms for its computations. In Section 3, the generalized SVD
(GSVD) and its extension to tensors (GTSVD) are outlined and it is shown how
to compute these decompositions. Two randomized algorithms are proposed in
Section 4 and the computer simulation results are reported in Section 6.
Finally, a conclusion is given in Section 7.

##  2 Basic definitions and concepts

We adopt the same notations used in [22] in this paper. So, to represent a
tensor, a matrix, and a vector we use an underlined bold capital letter, a
bold capital letter, and a bold lower letter. Slices are important subtensors
that are generated by fixing all but two modes. In particular, for a third-
order tensor $\underline{{\bf X}},$ the three types slices $\underline{{\bf
X}}(:,:,k),\,\underline{{\bf X}}(:,j,:),\,\underline{{\bf X}}(i,:,:)$ are
called frontal, lateral and horizontal slices. For convenience, sometimes in
the paper, we use the equivalent notations $\underline{{\bf
X}}_{i}\equiv\underline{{\bf X}}(;,:,i)$. Fibers are generated by fixing all
but one mode, so they are vectors. For a third-order tensor $\underline{{\bf
X}},$ the fiber $\underline{{\bf X}}(i,j,:)$ is called a tube. The notation
“conj” denotes the component-wise complex conjugate of a matrix. The Frobenius
norm of matrices or tensors is denoted by $\|.\|$. The notation $\|.\|_{2}$
stands for the Euclidean norm of vectors and the spectral norm of matrices.
For a positive definite matrix ${\bf S}$, the weighted inner product norm is
defined as $\|{\bf X}\|_{\bf S}=\sqrt{{\rm Tr}({\bf X}^{T}{\bf S}{\bf X})}$
where “${\rm Tr}$” is the trace operator. The mathematical expectation is
represented by $\mathbb{E}$. We now present the next definitions, which we
need in our presentation.

######  Definition 1.

(t-product) Let $\underline{\mathbf{X}}\in\mathbb{R}^{I_{1}\times I_{2}\times
I_{3}}$ and $\underline{\mathbf{Y}}\in\mathbb{R}^{I_{2}\times I_{4}\times
I_{3}}$, the tubal product (t-product)
$\underline{\mathbf{X}}*\underline{\mathbf{Y}}\in\mathbb{R}^{I_{1}\times
I_{4}\times I_{3}}$ is defined as follows

$\underline{\mathbf{C}}=\underline{\mathbf{X}}*\underline{\mathbf{Y}}={\rm
fold}\left({{\rm circ}\left(\underline{\mathbf{X}}\right){\rm
unfold}\left(\underline{\mathbf{Y}}\right)}\right),$ (1)

where

${\rm
circ}\left(\underline{\mathbf{X}}\right)=\begin{bmatrix}\underline{\mathbf{X}}(:,:,1)&\underline{\mathbf{X}}(:,:,I_{3})&\cdots&\underline{\mathbf{X}}(:,:,2)\\\
\underline{\mathbf{X}}(:,:,2)&\underline{\mathbf{X}}(:,:,1)&\cdots&\underline{\mathbf{X}}(:,:,3)\\\
\vdots&\vdots&\ddots&\vdots\\\
\underline{\mathbf{X}}(:,:,I_{3})&\underline{\mathbf{X}}(:,:,I_{3}-1)&\cdots&\underline{\mathbf{X}}(:,:,1)\end{bmatrix},$

and

${\rm
unfold}(\underline{\mathbf{Y}})=\begin{bmatrix}\underline{\mathbf{Y}}(:,:,1)\\\
\underline{\mathbf{Y}}(:,:,2)\\\ \vdots\\\
\underline{\mathbf{Y}}(:,:,I_{3})\end{bmatrix},\hskip
14.22636pt\underline{\mathbf{Y}}={\rm fold}\left({\rm
unfold}\left(\underline{\mathbf{Y}}\right)\right).$

######  Definition 2.

(Transpose) The transpose of a tensor
$\underline{\mathbf{X}}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is
denoted by $\underline{\mathbf{X}}^{T}\in\mathbb{R}^{I_{2}\times I_{1}\times
I_{3}}$ produced by applying the transpose to all frontal slices of the tensor
$\underline{\mathbf{X}}$ and reversing the order of the transposed frontal
slices from the second till the last one.

######  Definition 3.

(Identity tensor) Identity tensor
$\underline{\mathbf{I}}\in\mathbb{R}^{I_{1}\times I_{1}\times I_{3}}$ is a
tensor whose first frontal slice is an identity matrix of size $I_{1}\times
I_{1}$ and all other frontal slices are zero. It is easy to show
$\underline{{\bf I}}*\underline{{\bf X}}=\underline{{\bf X}}$ and
$\underline{{\bf X}}*\underline{{\bf I}}=\underline{{\bf X}}$ for all tensors
of conforming sizes.

######  Definition 4.

(Orthogonal tensor) We all that a tensor
$\underline{\mathbf{X}}\in\mathbb{R}^{I_{1}\times I_{1}\times I_{3}}$ is
orthogonal if
${\underline{\mathbf{X}}^{T}}*\underline{\mathbf{X}}=\underline{\mathbf{X}}*{\underline{\mathbf{X}}^{T}}=\underline{\mathbf{I}}$.

######  Definition 5.

(f-diagonal tensor) If all frontal slices of a tensor are diagonal then the
tensor is called an f-diagonal tensor.

######  Definition 6.

(Moore-Penrose pseudoinverse of a tensor) Let $\underline{\bf
X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ be given. The Moore-Penrose
(MP) pseudoinverse of the tensor $\underline{\bf X}$ is denoted by
$\underline{\bf X}^{{\dagger}}\in\mathbb{R}^{I_{2}\times I_{1}\times I_{3}}$
is a unique tensor satisfying the following four equations:

$\displaystyle\underline{\bf X}^{{\dagger}}*\underline{\bf X}*\underline{\bf
X}^{{\dagger}}=\underline{\bf X}^{{\dagger}},\quad\underline{\bf
X}*\underline{\bf X}^{{\dagger}}*\underline{\bf X}=\underline{\bf X},$
$\displaystyle(\underline{\bf X}*\underline{\bf
X}^{{\dagger}})^{T}=\underline{\bf X}*\underline{\bf
X}^{{\dagger}},\quad(\underline{\bf X}^{{\dagger}}*\underline{\bf
X})^{T}=\underline{\bf X}^{{\dagger}}*\underline{\bf X}.$

The MP pseudoinverse of a tensor can also be computed in the Fourier domain
and this is shown in Algorithm 2.

The inverse of a tensor is a special case of the MP pseudoinverse of tensors.
The inverse of $\underline{\bf X}\in\mathbb{R}^{I_{1}\times I_{1}\times
I_{3}}$ is denoted by $\underline{\bf X}^{-1}\in\mathbb{R}^{I_{1}\times
I_{1}\times I_{3}}$ is a unique tensor satisfying $\underline{\bf
X}*\underline{\bf X}^{-1}=\underline{\bf X}^{-1}*\underline{\bf
X}=\underline{\bf I},$ where $\underline{\bf I}\in\mathbb{R}^{I_{1}\times
I_{1}\times I_{3}}$ is the identity tensor The inverse of a tensor can be also
computed in the Fourier domain by replacing the MATLAB command “inv” in Line 3
of Algorithm 2.

Input : Two data tensors
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}},\,\,\underline{\mathbf{Y}}\in{\mathbb{R}^{{I_{2}}\times{I_{4}}\times{I_{3}}}}$

Output : t-product
$\underline{\mathbf{C}}=\underline{\mathbf{X}}*\underline{\mathbf{Y}}\in\mathbb{R}^{I_{1}\times
I_{4}\times I_{3}}$

1 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right)$;

2 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{Y}},[],3}\right)$;

3 for  _$i=1,2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil$_ do

4
$\widehat{\underline{\mathbf{C}}}\left({:,:,i}\right)=\widehat{\underline{\mathbf{X}}}\left({:,:,i}\right)\,\widehat{\underline{\mathbf{Y}}}\left({:,:,i}\right)$;

5

6 end for

7for  _$i=\lceil\frac{I_{3}+1}{2}\rceil+1,\ldots,I_{3}$_ do

8 $\widehat{\underline{\mathbf{C}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{C}}}\left({:,:,I_{3}-i+2}\right))$;

9 end for

$\underline{\mathbf{C}}={\rm
ifft}\left({\widehat{\underline{\mathbf{C}}},[],3}\right)$;

Algorithm 1 The t-product of two tensors [15, 23]

Input : The data tensor
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}}$

Output : Moore-Penrose pseudoinvers
$\underline{\mathbf{X}}^{{\dagger}}\in\mathbb{R}^{I_{2}\times I_{1}\times
I_{3}}$

1 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right)$;

2 for  _$i=1,2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil$_ do

3 $\widehat{\underline{\mathbf{C}}}\left({:,:,i}\right)={\rm
pinv}\,(\widehat{{\bf X}}(:,:,i))$;

4

5 end for

6for  _$i=\lceil\frac{I_{3}+1}{2}\rceil+1,\ldots,I_{3}$_ do

7 $\widehat{\underline{\mathbf{C}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{C}}}\left({:,:,I_{3}-i+2}\right))$;

8 end for

$\underline{\mathbf{X}}^{{\dagger}}={\rm
ifft}\left({\widehat{\underline{\mathbf{C}}},[],3}\right)$;

Algorithm 2 Fast Moore-Penrose pseudoinvers computation of the tensor
$\underline{\bf X}$

###  2.1 Tensor SVD (t-SVD), Tensor QR (t-QR) decomposition and Tensor LU
(t-LU) decomposition

The classical matrix decompositions such as QR, LU, and SVD can be
straightforwardly generalized to tenors based on the t-product. Given a tensor
$\underline{{\bf X}}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}},$ the tensor
QR (t-QR) decomposition represents the tensor $\underline{\bf X}$ as
$\underline{{\bf X}}=\underline{{\bf Q}}*\underline{{\bf R}}$ and can be
computed through Algorithm 3. By a slight modification of Algorithm 3, the
tensor LU (t-LU) decomposition and tensor SVD (t-SVD) can be computed. More
precisely, in line 4 of Algorithm 3, we replace the LU decomposition and the
SVD of frontal slices $\widehat{\underline{{\bf
X}}}(:,:,i),\,i=1,2,\ldots,I_{3},$ instead of the QR decomposition.

The tensor SVD (t-SVD) expresses a tensor as the t-product of three tensors.
The first and last tensors are orthogonal while the middle tensor if an
f-diagonal tensor. Let $\underline{{\bf X}}\in\mathbb{R}^{I_{1}\times
I_{2}\times I_{3}}$, then the t-SVD gives the following model:

$\underline{{\bf X}}\approx\underline{\bf U}*\underline{\bf S}*\underline{\bf
V}^{T},$

where $\underline{\bf U}\in\mathbb{R}^{I_{1}\times R\times
I_{3}},\,\underline{\bf S}\in\mathbb{R}^{R\times R\times I_{3}},$ and
$\underline{\bf V}\in\mathbb{R}^{I_{2}\times R\times I_{3}}$. The tensors
$\underline{\bf U}$ and $\underline{\bf V}$ are orthogonal, while the tensor
$\underline{\bf S}$ is f-diagonal. For the generalization of the t-SVD to
tensors of order higher than three is done in [24]. The t-SVD can be computed
via Algorithm 4.

Input : The data tensor
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}}$

Output : The t-QR decomposition of the tensor $\underline{\bf X}$

1 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right)$;

2 for  _$i=1,2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil$_ do

3 $[\widehat{\underline{\mathbf{Q}}}\left({:,:,i}\right),\widehat{{\bf
R}}(:,:,i)]={\rm qr}\,(\widehat{{\bf X}}(:,:,i),0)$;

4

5 end for

6for  _$i=\lceil\frac{I_{3}+1}{2}\rceil+1,\ldots,I_{3}$_ do

7 $\widehat{\underline{\mathbf{Q}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{Q}}}\left({:,:,I_{3}-i+2}\right))$;

8 $\widehat{\underline{\mathbf{R}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{R}}}\left({:,:,I_{3}-i+2}\right))$;

9 end for

10$\underline{\mathbf{Q}}={\rm
ifft}\left({\widehat{\underline{\mathbf{Q}}},[],3}\right)$;

$\underline{\mathbf{R}}={\rm
ifft}\left({\widehat{\underline{\mathbf{R}}},[],3}\right)$;

Algorithm 3 The t-QR decomposition of the tensor $\underline{\bf X}$

Input : The data tensor
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}}$
and a target tubal rank $R$

Output : The truncated t-SVD of the tensor $\underline{\bf X}$

1 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right)$;

2 for  _$i=1,2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil$_ do

3
$[\widehat{\underline{\mathbf{U}}}\left({:,:,i}\right),\widehat{\underline{\bf
S}}(:,:,i),\widehat{{\bf V}}(:,:,i)]={\rm svds}\,(\widehat{{\bf
X}}(:,:,i),R)$;

4

5 end for

6for  _$i=\lceil\frac{I_{3}+1}{2}\rceil+1,\ldots,I_{3}$_ do

7 $\widehat{\underline{\mathbf{U}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{U}}}\left({:,:,I_{3}-i+2}\right))$;

8
$\widehat{\underline{\mathbf{S}}}\left({:,:,i}\right)=\widehat{\underline{\mathbf{S}}}\left({:,:,I_{3}-i+2}\right)$;

9 $\widehat{\underline{\mathbf{V}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{V}}}\left({:,:,I_{3}-i+2}\right))$;

10 end for

$\underline{\mathbf{U}}_{R}={\rm
ifft}\left({\widehat{\underline{\mathbf{U}}},[],3}\right)$;
$\underline{\mathbf{S}}_{R}={\rm
ifft}\left({\widehat{\underline{\mathbf{S}}},[],3}\right)$;
$\underline{\mathbf{V}}_{R}={\rm
ifft}\left({\widehat{\underline{\mathbf{V}}},[],3}\right)$

Algorithm 4 The truncated t-SVD decomposition of the tensor $\underline{\bf
X}$

##  3 Generalized singular value decomposition (GSVD) and its extension to
tensors based on the t-product (GTSVD)

In this section, the GSVD and its extension to tensors based on the t-product
are introduced. The GSVD is a generalized version of the classical SVD, which
is applied to a pair of matrices. The SVD was generalized in [2] from two
different perspectives. More precisely, from the SVD, it is known that each
matrix ${\bf X}\in\mathbb{R}^{I_{1}\times I_{2}}$ can be decomposed in the
form ${\bf X}={\bf U}{\bf S}{\bf V}^{T}$ where ${\bf
U}\in\mathbb{R}^{I_{1}\times I_{1}}$ and ${\bf V}\in\mathbb{R}^{I_{2}\times
I_{2}}$ are orthogonal matrices and the ${\bf S}={\rm
diag}(\sigma_{1},\sigma_{2},\ldots,\sigma_{\min\\{I_{1},I_{2}\\}})\in\mathbb{R}^{I_{1}\times
I_{2}}$ is a diagonal matrix with singular values
$\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{R}>\sigma_{R+1}=\ldots=\sigma_{\min\\{I_{1},I_{2}\\}}=0$
and ${\rm rank}({\bf X})=R.$ Denoting the set of singular values of the matrix
${\bf X}$ as $\sigma({\bf
X})=\\{\sigma_{1},\sigma_{2},\ldots,\sigma_{\min\\{I_{1},I_{2}\\}}\\}$, it is
known that

$\displaystyle\sigma_{i}\in\sigma({\bf X})\longrightarrow{\rm det}({\bf
X}^{T}{\bf X}-\sigma_{i}^{2}{\bf I})=0,$ (2)
$\displaystyle\sigma_{i}\in\sigma({\bf
X})\longrightarrow\sigma_{i}\textrm{\,is a stationary point of\,}\frac{\|{\bf
X}{\bf z}\|_{2}}{\|{\bf z}\|_{2}}.$ (3)

Based on (2) and 3, the SVD was generalized in the following straightforward
ways:

$\displaystyle\textrm{Find\,}\sigma\geq 0\,\textrm{such that\,}{\rm det}({\bf
X}^{T}{\bf X}-\sigma^{2}{\bf Y}^{T}{\bf Y})=0,$ (4) $\displaystyle\textrm{Find
the stationary values of\,}\frac{\|{\bf X}{\bf z}\|_{\bf S}}{\|{\bf z}\|_{\bf
T}},$ (5)

where ${\bf Y}\in\mathbb{R}^{I_{3}\times I_{2}}$ is an arbitrary matrix and
${\bf S}\in\mathbb{R}^{I_{2}\times I_{2}}$ and ${\bf
T}\in\mathbb{R}^{I_{1}\times I_{1}}$ are positive definitive matrices. Let us
denote the set of all points satisfying (4) as $\sigma({\bf X},{\bf
Y})=\\{\sigma|\sigma\geq 0,\,{\rm det}({\bf X}^{T}{\bf X}-\sigma^{2}{\bf
Y}^{T}{\bf Y})=0\\},$ which are called ${\bf Y}$-singular values of the matrix
${\bf X}$. It was shown in [2] that given ${\bf X}\in\mathbb{R}^{I_{1}\times
I_{2}},\,\,I_{1}\geq I_{2}$ and ${\bf Y}\in\mathbb{R}^{I_{3}\times
I_{2}},\,I_{3}\geq I_{2}$ there exist orthogonal matrices ${\bf
U}\in\mathbb{R}^{I_{1}\times I_{1}},\,{\bf V}\in\mathbb{R}^{I_{3}\times
I_{3}}$ and a nonsingular matrix ${\bf Z}\in\mathbb{R}^{I_{2}\times I_{2}}$
such that

$\displaystyle{\bf U}^{T}{\bf X}{\bf Z}$ $\displaystyle=$ $\displaystyle{\rm
diag}(\alpha_{1},\cdots,\alpha_{I_{2}}),\quad\alpha_{i}\in[0,1]$ (6)
$\displaystyle{\bf V}^{T}{\bf Y}{\bf Z}$ $\displaystyle=$ $\displaystyle{\rm
diag}(\beta_{1},\cdots,\beta_{I_{2}}),\quad\beta_{i}\in[0,1]$ (7)

where $\alpha_{i}^{2}+\beta_{i}^{2}=1$ with the ratios $\beta_{i}/\alpha_{i}$
of nondecreasing order for $i=1,2,\ldots,I_{2}$. A more general SVD was
proposed in [3] where a more computationally stable algorithm was also
developed to compute it. In the following, the latter GSVD is introduced,
which will be considered in our paper.

######  Theorem 1.

[3] Let two matrices ${\bf X}\in\mathbb{R}^{I_{1}\times I_{2}}$ and ${\bf
Y}\in\mathbb{R}^{I_{3}\times I_{2}}$ be given and assume that the SVD of the
matrix ${\bf C}=\begin{bmatrix}{\bf X}\\\ {\bf Y}\end{bmatrix}$ is

$\displaystyle{\bf E}^{T}{\bf C}{\bf Z}=\begin{bmatrix}{\bf\Gamma}&{\bf 0}\\\
{\bf 0}&{\bf 0}\end{bmatrix},$ (8)

with the unitary matrices ${\bf
E}\in\mathbb{R}^{(I_{1}+I_{3})\times(I_{1}+I_{3})},\,{\bf
Z}\in\mathbb{R}^{I_{2}\times I_{2}}$ and a diagonal matrix
${\bf\Gamma}\in\mathbb{R}^{k\times k}$. Then, there exist matrices ${\bf
U}\in\mathbb{R}^{I_{1}\times I_{1}},\,{\bf V}\in\mathbb{R}^{I_{3}\times
I_{3}}$ and ${\bf W}\in\mathbb{R}^{k\times k}$ such that

$\displaystyle{\bf U}^{T}{\bf X}{\bf Z}={\bf\Sigma}_{\bf X}({\bf
W}^{T}{\bf\Gamma},{\bf 0}),\quad{\bf V}^{T}{\bf Y}{\bf Z}={\bf\Sigma}_{\bf
Y}({\bf W}^{T}{\bf\Gamma},{\bf 0}),$ (9)

where ${\bf\Sigma}_{\bf X}\in\mathbb{R}^{I_{1}\times k}$ and ${\bf\Sigma}_{\bf
Y}\in\mathbb{R}^{I_{3}\times k}$ are defined as follows:

$\displaystyle{\bf\Sigma}_{X}=\begin{bmatrix}{\bf I}_{\bf X}&&\\\ &{\bf
S}_{\bf X}&\\\ &&{\bf 0}_{\bf
X}\end{bmatrix},\quad{\bf\Sigma}_{Y}=\begin{bmatrix}{\bf 0}_{\bf Y}&&\\\ &{\bf
S}_{\bf Y}&\\\ &&{\bf I}_{\bf Y}\end{bmatrix}.$ (10)

Here, ${\bf I}_{\bf X}\in\mathbb{R}^{c\times c}$ and ${\bf I}_{\bf
Y}\in\mathbb{R}^{(k-c-d)\times(k-c-d)}$ are identity matrices, ${\bf 0}_{\bf
X}\in\mathbb{R}^{(I_{1}-c-d)\times(k-c-d)}$ and ${\bf 0}_{\bf
Y}\in\mathbb{R}^{(I_{3}-k+c)\times c}$ are zero matrices that may have no
columns/rows and ${\bf S}_{\bf X}\in\mathbb{R}^{d\times d},\,{\bf S}_{\bf
Y}\in\mathbb{R}^{d\times d}$ are diagonal matrices with diagonal elements
$1>\alpha_{c+1}\geq\cdots\geq\alpha_{c+d}>0$ and
$1<\beta_{c+1}\leq\cdots\leq\beta_{c+d}>0$, respectively and
$\alpha_{i}^{2}+\beta_{i}^{2}=1$ for $c+1\leq i\leq c+d$. Note that $c$ and
$d$ are internally defined by the matrices ${\bf X}$ and ${\bf Y}$.

It is not difficult to check that (9) is reduced to

$\displaystyle{\bf U}^{T}{\bf X}{\bf R}^{-1}=({\bf\Sigma}_{\bf X},{\bf
0}),\quad{\bf V}^{T}{\bf Y}{\bf R}^{-1}=({\bf\Sigma}_{\bf Y},{\bf 0}),$ (11)

for ${\bf R}^{-1}$ defined as follows

${\bf R}^{-1}={\bf Z}\begin{bmatrix}{\bf\Gamma}^{-1}{\bf W}&{\bf 0}\\\ {\bf
0}&{\bf I}_{I_{2}-k}\end{bmatrix},$

and if the matrix ${\bf C}$ is of full rank, then the zero blocks on the
right-hand sides of (11) are removed. As we see, the first formulation (4) of
the GSVD deals with two matrices ${\bf X}$ and ${\bf Y},$ and provides a
decomposition of the form (11) with the same right-hand side matrix ${\bf
R}^{-1}$. There is a new GSVD corresponding to (5) given two positive definite
matrices ${\bf S}$ and ${\bf T}$, see [2] for details. In this paper, we only
consider the GSVD formulated in (2) and its extension to the tensor case based
on the t-product.

The GSVD can be analogously extended to tensors based on the t-product [20,
21]. Let $\underline{\bf X},\,\underline{\bf Y}$ be two given tensors with the
same number of lateral slices. Then the Generalized tensor SVD (GTSVD)
decomposes the tensors $\underline{\bf X}\in\mathbb{R}^{I_{1}\times
I_{2}\times I_{3}},\,\,I_{1}\geq I_{2},\,\underline{\bf
Y}\in\mathbb{R}^{I_{4}\times I_{2}\times I_{3}},\,I_{4}\geq I_{2},$ jointly in
the following form:

$\displaystyle\underline{\bf X}$ $\displaystyle=$ $\displaystyle\underline{\bf
U}*\underline{\bf C}*\underline{\bf Z},$ (12) $\displaystyle\underline{\bf Y}$
$\displaystyle=$ $\displaystyle\underline{\bf V}*\underline{\bf
S}*\underline{\bf Z},$ (13)

where $\underline{\bf U}\in\mathbb{R}^{I_{1}\times I_{1}\times
I_{3}},\,\underline{\bf V}\in\mathbb{R}^{I_{4}\times I_{4}\times
I_{3}},\,\underline{\bf C}\in\mathbb{R}^{I_{1}\times I_{2}\times
I_{3}},\,\underline{\bf S}\in\mathbb{R}^{I_{4}\times I_{2}\times
I_{3}},\,\underline{\bf Z}\in\mathbb{R}^{I_{2}\times I_{2}\times I_{3}}$. Note
that the tensors $\underline{\bf C}$ and $\underline{\bf S}$ are f-diagonal
and the tensors $\underline{\bf U}$ and $\underline{\bf V}$ are orthogonal and
$\underline{\bf Z}$ is nonsingular. The procedure of the computation of the
GTSVD is presented in Algorithm 5. We need to apply the classical GSVD (lines
3-5) to the first $\lceil\frac{I_{3}+1}{2}\rceil$ frontal slices of the
tensors $\underline{\bf X}$ and $\underline{\bf Y}$ in the Fourier domain and
the rest of the slices are computed easily (Lines 6-12).

The computation of the GSVD or the GTSVD for large-scale matrices/tensors
involves the computation of the SVD of some large matrices. So, it is
computationally demanding and requires huge memory and resources. In recent
years, the idea of randomization has been utilized to accelerate the
computation of the GSVD. In [8], a randomized algorithm is proposed for the
GSVD of the form (4) while [9] proposes a randomized algorithm for the GSVD of
form 5. The key idea is to employ the random projection method for fast
computation of the SVD, which is required in the process of computing the
GSVD.

Input : The data tensors
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}}$
and
$\underline{\mathbf{Y}}\in{\mathbb{R}^{{I_{4}}\times{I_{2}}\times{I_{3}}}}$

Output : The generalized t-SVD (GTSVD) of $\underline{\mathbf{X}}$ and
$\underline{\bf Y}$ as
$\underline{\mathbf{X}}=\underline{\mathbf{U}}*\underline{\mathbf{C}}*\underline{\mathbf{Z}}$
and
$\underline{\mathbf{Y}}=\underline{\mathbf{V}}*\underline{\mathbf{S}}*\underline{\mathbf{Z}}$

1 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right)$;

2 $\widehat{\underline{\mathbf{Y}}}={\rm
fft}\left({\underline{\mathbf{Y}},[],3}\right)$;

3 for  _$i=1,2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil$_ do

4
$[\widehat{\underline{\mathbf{U}}}_{i},\,\widehat{\underline{\mathbf{V}}}_{i},\,\widehat{\underline{\mathbf{Z}}}_{i},\,\widehat{\underline{\mathbf{C}}}_{i},\,\widehat{\underline{\mathbf{S}}}_{i}]={\rm
GSVD}\,(\widehat{\underline{{\bf X}}}(:,:,i),\widehat{\underline{{\bf
Y}}}(:,:,i))$;

5

6 end for

7for  _$i=\lceil\frac{I_{3}+1}{2}\rceil+1,\ldots,I_{3}$_ do

8 $\widehat{\underline{\mathbf{U}}}_{i}={\rm
conj}(\widehat{\underline{\mathbf{U}}}_{I_{3}-i+2})$;

9 $\widehat{\underline{\mathbf{V}}}_{i}={\rm
conj}({\underline{\mathbf{V}}}_{I_{3}-i+2})$;

10 $\widehat{\underline{\mathbf{Z}}}_{i}={\rm
conj}(\widehat{\underline{\mathbf{Z}}}_{I_{3}-i+2})$;

11
$\widehat{\underline{\mathbf{C}}}_{i}=\widehat{\underline{\mathbf{C}}}_{I_{3}-i+2}$;

12
$\widehat{\underline{\mathbf{S}}}_{i}=\widehat{\underline{\mathbf{S}}}_{I_{3}-i+2}$;

13 end for

${\underline{\mathbf{U}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{U}}},[],3}\right)$;
$\,{\underline{\mathbf{V}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{V}}},[],3}\right)$;
$\,{\underline{\mathbf{Z}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{Z}}},[],3}\right)$;
$\,{\underline{\mathbf{C}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{C}}},[],3}\right)$;
$\,{\underline{\mathbf{S}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{S}}},[],3}\right)$;

Algorithm 5 Generalized t-SVD of $\underline{\bf X}$ and $\underline{\bf Y}$

##  4 Proposed fast randomized algorithms for computation of the GTSVD

In this section, we propose two randomized variants of the GTSVD Algorithm 5.
The first proposed randomized algorithm is a naive modification of Algorithm 5
where we can replace the deterministic GSVD with the randomized counterpart
developed in [8]. This idea is presented in Algorithm 6. Here, we can use the
oversampling and the power iteration methods to improve the accuracy of the
singular values of the frontal slices do not decay sufficiently fast [7]. The
second proposed randomized algorithm is presented in Algorithm 7 where we
first make a reduction on the given data tensors by multiplying (t-product)
them with random tensors to capture their important actions. Similar to
Algorithm 6, if the frontal slices of a given tensor do not have a fast decay,
the power iteration technique and the oversampling method should be used to
better capture their ranges. More precisely, the random projection stages in
Algorithm 7 (Lines 1-2) are replaced with the following computations

$\displaystyle\underline{\bf W}_{1}=(\underline{\bf X}*\underline{\bf
X}^{T})^{q}*\underline{\bf X}*\underline{\bf\Omega}_{1},$ (14)
$\displaystyle\underline{\bf W}_{2}=(\underline{\bf Y}*\underline{\bf
Y}^{T})^{q}*\underline{\bf Y}*\underline{\bf\Omega}_{2}.$ (15)

In practice, for the computations (14)-(15) to be stable, we employ the t-QR
decomposition or the t-LU decomposition or a combination of them [7, 25, 13].
Then, by applying the t-QR algorithm to the mentioned compressed tensors, we
can obtain orthonormal bases for them, which are used to get small sketch
tensors by projecting the original data tensors onto the compressed tensors.
Since the sizes of the sketch tensors are smaller than the original ones, the
deterministic algorithms can be used to compute their GTSVDs. Finally, the
GTSVD of the original data tensor can be recovered from the GTSVD of the
compressed tensors. Note that in Algorithms 6 and 7, we need the tubal rank as
input, however, this can be numerically estimated for a given approximation
error bound. For example, we can use the randomized fixed-precision developed
in [26] and for the case of tensors, the randomized rank-revealing algorithm
proposed in [27, 28, 25] is applicable.

Input : The data tensors
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}}$
and
$\underline{\mathbf{Y}}\in{\mathbb{R}^{{I_{4}}\times{I_{2}}\times{I_{3}}}}$,
standard Gaussian matrices with oversampling parameter $p_{1}$ and $p_{2}$

Output : The generalized t-SVD of $\underline{\mathbf{X}}$ and $\underline{\bf
Y}$ as
$\underline{\mathbf{X}}=\underline{\mathbf{U}}*\underline{\mathbf{C}}*\underline{\mathbf{Z}}$
and
$\underline{\mathbf{Y}}=\underline{\mathbf{V}}*\underline{\mathbf{S}}*\underline{\mathbf{Z}}$

1 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right)$;

2 $\widehat{\underline{\mathbf{Y}}}={\rm
fft}\left({\underline{\mathbf{Y}},[],3}\right)$;

3 for  _$i=1,2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil$_ do

4
$[\widehat{\underline{\mathbf{U}}}_{i},\,\widehat{\underline{\mathbf{V}}}_{i},\,\widehat{\underline{\mathbf{Z}}}_{i},\,\widehat{\underline{\mathbf{C}}}_{i},\,\widehat{\underline{\mathbf{S}}}_{i}]={\rm
Randomized\,\,GSVD}\,(\widehat{\underline{{\bf
X}}}(:,:,i),\widehat{\underline{{\bf Y}}}(:,:,i))$;

5

6 end for

7for  _$i=\lceil\frac{I_{3}+1}{2}\rceil+1,\ldots,I_{3}$_ do

8 $\widehat{\underline{\mathbf{U}}}_{i}={\rm
conj}(\widehat{\underline{\mathbf{U}}}_{I_{3}-i+2})$;

9 $\widehat{\underline{\mathbf{V}}}_{i}={\rm
conj}({\underline{\mathbf{V}}}_{I_{3}-i+2})$;

10 $\widehat{\underline{\mathbf{Z}}}_{i}={\rm
conj}(\widehat{\underline{\mathbf{Z}}}_{I_{3}-i+2})$;

11
$\widehat{\underline{\mathbf{C}}}_{i}=\widehat{\underline{\mathbf{C}}}_{I_{3}-i+2}$;

12
$\widehat{\underline{\mathbf{S}}}_{i}=\widehat{\underline{\mathbf{S}}}_{I_{3}-i+2}$;

13 end for

$\widehat{\underline{\mathbf{U}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{U}}},[],3}\right)$;
$\,\widehat{\underline{\mathbf{V}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{V}}},[],3}\right)$;
$\,\widehat{\underline{\mathbf{Z}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{Z}}},[],3}\right)$;
$\,\widehat{\underline{\mathbf{C}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{C}}},[],3}\right)$;
$\,\widehat{\underline{\mathbf{S}}}={\rm
ifft}\left({\widehat{\underline{\mathbf{S}}},[],3}\right)$;

Algorithm 6 The proposed randomized GTSVD I

Input : The data tensors
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}}$
and
$\underline{\mathbf{Y}}\in{\mathbb{R}^{{I_{4}}\times{I_{2}}\times{I_{3}}}}$,
standard Gaussian tensors
$\underline{\bf\Omega}_{1},\,\underline{\bf\Omega}_{2}$ with corresponding
oversampling parameters $p_{1}$ and $p_{2}$

Output : The generalized t-SVD of $\underline{\mathbf{X}}$ and $\underline{\bf
Y}$ as
$\underline{\mathbf{X}}=\underline{\mathbf{U}}*\underline{\mathbf{C}}*\underline{\mathbf{Z}}$
and
$\underline{\mathbf{Y}}=\underline{\mathbf{V}}*\underline{\mathbf{S}}*\underline{\mathbf{Z}}$

1 $\underline{\bf W}_{1}=\underline{\bf X}*\underline{\bf\Omega}_{1}$;

2 $\underline{\bf W}_{2}=\underline{\bf Y}*\underline{\bf\Omega}_{2}$;

3 $[\underline{\bf Q}_{1},\sim]=t$-${\rm QR}(\underline{\bf W}_{1})$;

4 $[\underline{\bf Q}_{2},\sim]=t$-${\rm QR}(\underline{\bf W}_{2})$;

5 $[\widehat{\underline{\bf U}},\underline{\bf C},\underline{\bf Z}]={\rm
GTSVD}(\underline{\bf Q}_{1}^{T}*\underline{\bf X})$;

6 $[\widehat{\underline{\bf V}},\underline{\bf S},\sim]={\rm
GTSVD}(\underline{\bf Q}_{2}^{T}*\underline{\bf Y})$;

7 ${\underline{\bf U}}={\underline{\bf Q}}_{1}*\widehat{\underline{\bf U}}$;

${\underline{\bf V}}={\underline{\bf Q}}_{2}*\widehat{\underline{\bf V}}$;

Algorithm 7 The proposed randomized GTSVD II

##  5 Error Analysis

We provide the average/expected error bounds of the approximations obtained by
the proposed randomized algorithms in this section. Let us first partition the
GSVD in (11) where ${\bf Z}={\bf R}^{-1}$ in the following form

$\displaystyle{\bf X}={\bf U}\begin{bmatrix}{\bf\Sigma}_{{\bf X}_{1}}&{\bf
0}&{\bf 0}\\\ {\bf 0}&{\bf\Sigma}_{{\bf X}_{2}}&{\bf 0}\\\
\end{bmatrix}\begin{bmatrix}{\bf Z}_{1}\\\ {\bf Z}_{2}\\\ {\bf
Z}_{3}\end{bmatrix},\quad{\bf Y}={\bf V}\begin{bmatrix}{\bf\Sigma}_{{\bf
Y}_{1}}&{\bf 0}&{\bf 0}\\\ {\bf 0}&{\bf\Sigma}_{{\bf Y}_{2}}&{\bf 0}\\\
\end{bmatrix}\begin{bmatrix}\widehat{{\bf Z}}_{1}\\\ \widehat{{\bf Z}}_{2}\\\
{\bf Z}_{3}\end{bmatrix}$ (16)

where ${\bf\Sigma}_{{\bf X}_{1}}\in\mathbb{R}^{r_{1}\times
r_{1}},\,{\bf\Sigma}_{{\bf
X}_{2}}\in\mathbb{R}^{(I_{1}-r_{1})\times(k-r_{1})},\,{\bf\Sigma}_{{\bf
Y}_{1}}\in\mathbb{R}^{(I_{3}-r_{2})\times(k-r_{1})},\,{\bf\Sigma}_{{\bf
Y}_{2}}\in\mathbb{R}^{r_{2}\times r_{2}},\,{\bf
Z}_{1}\in\mathbb{R}^{r_{1}\times I_{2}},\,{\bf
Z}_{2}\in\mathbb{R}^{(k-r_{1})\times I_{2}},\,\widehat{{\bf
Z}}_{1}\in\mathbb{R}^{(k-r_{2})\times I_{3}},\,\widehat{{\bf
Z}}_{2}\in\mathbb{R}^{r_{2}\times I_{2}}$ and ${\bf
Z}_{3}\in\mathbb{R}^{(I_{2}-k)\times I_{2}}$. Also, consider the standard
Gaussian matrices ${\bf\Phi}\in\mathbb{R}^{I_{2}\times(r_{1}+p_{1})}$ and
${\bf\Phi}\in\mathbb{R}^{I_{2}\times(r_{2}+p_{2})}$ for the given target ranks
$r_{1},\,r_{2}$ and the oversampling parameters $p_{1},\,p_{2}$. We first
start with the following theorem that gives the average error bound of an
approximation yielded by the random projection method for the computation of
the GSVD [8].

######  Theorem 2.

[8] Given two matrices ${\bf X}\in\mathbb{R}^{I_{1}\times I_{2}},\,{\bf
Y}\in\mathbb{R}^{I_{3}\times I_{2}}$ and consider the target ranks
$r_{1},\,r_{2}$ and the oversampling parameters $p_{1},\,p_{2}$
($r_{1}+p_{1}\leq\min({I_{1},I_{2}}),\,r_{2}+p_{2}\leq\min({I_{3},I_{2}}))$.
Assume that two standard random matrices
${\bf\Phi}\in\mathbb{R}^{I_{2}\times(r_{1}+p_{1})}$ and
${\bf\Phi}\in\mathbb{R}^{I_{2}\times(r_{2}+p_{2})}$ are used to capture the
range of the matrices ${\bf X},\,{\bf Y}$ using ${\bf M}={\bf X}\Phi,\,{\bf
N}={\bf Y}\Psi,$ then

$\displaystyle\mathbb{E}(\|{\bf I}-{\bf M}{\bf M}^{{\dagger}}{\bf
X}\|)\leq{\eta^{\bf X}_{1}}{\alpha}_{r_{1}+1}+\eta^{{\bf
X}}_{2}\sqrt{\sum_{j>r_{1}}\alpha_{j}^{2}},$ (17)
$\displaystyle\mathbb{E}(\|{\bf I}-{\bf N}{\bf N}^{{\dagger}}{\bf
Y}\|)\leq\eta^{{\bf Y}}_{1}{\beta}_{k-r_{2}}+\eta^{{\bf
Y}}_{2}\sqrt{\sum_{j\leq k-r_{2}}\beta_{j}^{2}},$ (18)

where

$\displaystyle\eta^{\bf X}_{1}=\|{\bf Z}\|\left(1+\frac{\sigma_{1}({\bf
Z}_{2})}{\sigma_{r_{1}}({\bf
Z}_{1})}\right)+\sqrt{\frac{r_{1}}{p_{1}-1}\sum_{j=1}^{r_{1}}\frac{\sigma_{1}^{2}({\bf
Z}_{2})}{\sigma_{j}^{2}({\bf Z}_{1})}},\,\,\eta^{\bf X}_{2}=\|{\bf
Z}\|\frac{\|{\bf Z}_{2}\|_{F}}{\sigma_{r_{1}}({\bf
Z}_{1})}\frac{e\sqrt{r_{1}+p_{1}}}{p_{1}},$ (19) $\displaystyle\eta^{\bf
Y}_{1}=\|{\bf Z}\|\left(1+\frac{\sigma_{1}(\widehat{{\bf
Z}}_{1})}{\sigma_{r_{2}}(\widehat{{\bf
Z}}_{2})}\right)+\sqrt{\frac{r_{2}}{p_{2}-1}\sum_{j=1}^{r_{2}}\frac{\sigma_{1}^{2}(\widehat{{\bf
Z}}_{1})}{\sigma_{j}^{2}(\widehat{{\bf Z}}_{2})}},\,\,\eta^{\bf Y}_{2}=\|{\bf
Z}\|\frac{\|\widehat{{\bf Z}}_{1}\|_{F}}{\sigma_{r_{2}}(\widehat{{\bf
Z}}_{2})}\frac{e\sqrt{r_{2}+p_{2}}}{p_{2}}.$ (20)

The average error bound of the approximation obtained by Algorithms 6-7 are
provided in Theorem 3.

######  Theorem 3.

Let
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}}$
and
$\underline{\mathbf{Y}}\in{\mathbb{R}^{{I_{4}}\times{I_{2}}\times{I_{3}}}}$ be
given data tensors. Assume we use the standard random tensors
$\underline{\bf\Omega}\in\mathbb{R}^{I_{2}\times(r_{1}+p_{1})\times
I_{3}},\,\underline{\bf\Omega}_{2}\in\mathbb{R}^{I_{2}\times(r_{2}+p_{2})\times
I_{3}}$ for the reduction stage and let their compressed tensors are

$\underline{\bf W}_{1}=\underline{\bf
X}*\underline{\bf\Omega}_{1},\quad\underline{\bf W}_{2}=\underline{\bf
Y}*\underline{\bf\Omega}_{2},$

where $r_{1},\,r_{2}$ are the target tubal ranks and oversampling parameters
$p_{1},\,p_{2}$ are the oversampling parameters. Then, the GTSVDs computed by
Algorithms 6-7 provide the solutions with the following accuracies

$\displaystyle\mathbb{E}(\|({\bf I}-{\underline{\bf W}_{1}}*{\underline{\bf
W}_{1}^{{\dagger}}}*\underline{{\bf
X}})\|)\leq\left(\frac{1}{I_{3}}\sum_{i=1}^{I_{3}}\eta_{1}^{\widehat{\bf
X}^{(i)}}{\alpha}^{i}_{r_{1}+1}+\eta_{2}^{\widehat{\bf
X}^{(i)}}\sqrt{\sum_{j>r_{1}}({\alpha^{i}_{j}})^{2}}\right)^{1/2},$ (21)
$\displaystyle\mathbb{E}(\|({\bf I}-{\underline{\bf W}_{2}}*{\underline{\bf
W}_{2}^{{\dagger}}}*\underline{{\bf
Y}})\|)\leq\left(\frac{1}{I_{3}}\sum_{i=1}^{I_{3}}\eta_{1}^{\widehat{\bf
Y}^{(i)}}+\eta_{2}^{\widehat{\bf Y}^{(i)}}\sqrt{\sum_{j\leq
k-r_{2}}({\beta^{i}_{j}})^{2}}\right)^{1/2},$ (22)

where, according to (16), the GSVDs of the frontal slices $\widehat{\bf
X}^{(i)}=\widehat{{\bf X}}(:,:,i)$ and $\widehat{\bf Y}^{(i)}=\widehat{{\bf
Y}}(:,:,i),\,\,i=1,2,\ldots,I_{3}$ are partitioned as

$\displaystyle\widehat{\bf X}^{(i)}={\bf
U}^{(i)}\begin{bmatrix}{\bf\Sigma}_{{\bf X}^{(i)}_{1}}&{\bf 0}&{\bf 0}\\\ {\bf
0}&{\bf\Sigma}_{{\bf X}^{(i)}_{2}}&{\bf 0}\\\ \end{bmatrix}\begin{bmatrix}{\bf
Z}^{{(i)}}_{1}\\\ {\bf Z}^{{(i)}}_{2}\\\ {\bf
Z}^{{(i)}}_{3}\end{bmatrix},\,\,\widehat{\bf Y}^{(i)}={\bf
V}^{(i)}\begin{bmatrix}{\bf\Sigma}_{{\bf Y}^{{(i)}}_{1}}&{\bf 0}&{\bf 0}\\\
{\bf 0}&{\bf\Sigma}_{{\bf Y}^{{(i)}}_{2}}&{\bf 0}\\\
\end{bmatrix}\begin{bmatrix}\widehat{{\bf Z}}^{{(i)}}_{1}\\\ \widehat{{\bf
Z}}^{{(i)}}_{2}\\\ {\bf Z}^{{(i)}}_{3}\end{bmatrix},$ (23)

where $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right),\,\,\widehat{\underline{\mathbf{Y}}}={\rm
fft}\left({\underline{\mathbf{Y}},[],3}\right)$. Here, the quantities
$\eta_{1}^{\widehat{\bf X}^{(i)}},\,\eta_{2}^{\widehat{\bf X}^{(i)}}$ and
$\eta_{1}^{\widehat{\bf Y}^{(i)}},\,\eta_{2}^{\widehat{\bf Y}^{(i)}}$ are
defined analogously based on the matrices $\widehat{\bf
Z}_{1}^{(i)},\,\widehat{\bf Z}_{2}^{(i)},\,\widehat{\bf
Z}_{3}^{(i)},\,{\widehat{\bf Z}}_{1}^{(i)},$ and $\widehat{\bf Z}_{2}^{(i)},$
(replacing them in (19)-(20) instead of $\widehat{\bf Z}_{1},\,\widehat{\bf
Z}_{2},\,\widehat{\bf Z}_{3},\,{\widehat{\bf Z}}_{1},\,\widehat{\bf Z}_{2}$).
Also $\alpha_{j}^{i},\,\beta_{j}^{i},$
$i=1,2,\ldots,I_{3},\,j=1,2,\ldots,I_{2}$ are the elements of the diagonal
middle matrices obtained from the GSVD of the matrices $(\widehat{\bf
X}^{(i)},\widehat{\bf Y}^{(i)})$.

###### Proof.

We prove the theorem only for (21), and part (22) can be similarly proved.
Considering the linearity of the expectation operator and formula
(LABEL:eq_fou) we have

$\displaystyle\mathbb{E}\,(\|\underline{{\bf X}}-{{\bf{\bf
W}}_{1}}*{{\bf\underline{{\bf W}}}^{{\dagger}}_{1}}*\underline{{\bf
X}}\|^{2})\leq\hskip 56.9055pt$
$\displaystyle\frac{1}{I_{3}}\left(\sum_{i=1}^{I_{3}}\mathbb{E}\,\|\widehat{\bf
X}^{(i)}-\widehat{\bf W}_{1}^{(i)}\widehat{\bf
W}_{1}^{(i)\,{\dagger}}\widehat{\bf X}^{(i)}\|^{2}\right),$ (24)

where $\widehat{\bf X}^{(i)}=\widehat{{\bf X}}(:,:,i)$ and $\widehat{\bf
W}_{1}^{(i)}=\widehat{{\bf W}}_{1}(:,:,i)$. Now, if we use Theorem 2 to bound
each term of (5) and use the Holder’s identity

$\mathbb{E}({\,\|\underline{\bf X}-\underline{{\bf W}}_{1}*\underline{{\bf
W}}_{1}^{{\dagger}}*\underline{\bf
X}\|})\leq(\left(\mathbb{E}({\,\|\underline{\bf X}-\underline{\bf
W}_{1}*\underline{\bf W}_{1}^{{\dagger}}*\underline{\bf
X}\|^{2}})\right)^{1/2},$

the proof is concluded. ∎

Figure 1: Example 1. Running time comparison of different algorithms.

##  6 Experimental Results

In this section, we conduct several simulations to show the efficiency of the
proposed algorithms and their superiority over the baseline algorithm. We have
used Matlab and some functions of the toolbox:
https://github.com/canyilu/Tensor-tensor-product-toolbox to implement the
proposed algorithm using a laptop computer with 2.60 GHz Intel(R) Core(TM)
i7-5600U processor and 8GB memory. The algorithms are compared in terms of
running time and relative error defined as follows

$\frac{\|\underline{\mathbf{X}}-\underline{\mathbf{U}}*\underline{\mathbf{C}}*\underline{\mathbf{Z}}\|_{F}+\|\underline{\mathbf{Y}}-\underline{\mathbf{V}}*\underline{\mathbf{S}}*\underline{\mathbf{Z}}\|_{F}}{\|\underline{\bf
X}\|_{F}+\|\underline{\bf Y}\|_{F}}.$

######  Example 1.

Let us generate random data tensors $\underline{\bf X}$ and $\underline{\bf
Y}$ with zero mean and unit variance of size $n\times n\times n$ with the
tubal rank 50, where $n=200,300,400,500$. Then, the basic TGSVD and two
proposed randomized TGSVD algorithms (Algorithms 6 and 6) are applied to the
mentioned data tensor. We set the oversampling as $50$ in both Algorithms 6
and 7. The running times of the algorithms are shown in Figure 1 and the
corresponding relative errors achieved by them are reported in Table 1. From
Figure 1, for $n=300,400,500$ we achieve $\times 55,\,\times 37.9,$ and
$\times 11.01$ speed-up, respectively. So, in all scenarios, we have more than
one order of magnitude acceleration. Also from Table 1, we see that the
difference between the relative errors of the algorithms is negligible. So, we
can provide satisfying results, in much less time than the baseline algorithm
5. This shows the superiority of the proposed algorithms compared to the
baseline methods for handling large-scale data tensors.

Table 1: Relative error comparison of different algorithms. | $n=300$ | $n=400$ | $n=500$
---|---|---|---
Algorithm 5 | 9.8833e-19 | 6.7230e-19 | 7.4645e-19
Proposed Algorithm 6 | 5.4381e-18 | 2.4254e-18 | 2.4809e-18
Proposed Algorithm 7 | 9.8108e-18 | 9.5354e-18 | 2.2562e-18

######  Example 2.

In this example, we consider the following data tensors

  * 1.

$\underline{\bf X}(i,j,k)=\frac{1}{\sqrt{i^{2}+j^{2}+k^{2}}}$

  * 2.

$\underline{\bf Y}(i,j,k)=\frac{1}{(i^{3}+j^{3}+k^{3})^{1/3}}$

It is easy to check that these data tensors have low tubal-rank structures. We
set the oversampling parameter to $50$ and the tubal rank $R=50$ in Algorithms
6 and 7 and apply them to the mentioned data tensors. The execution times of
the proposed algorithms and the baseline Algorithm 5 are reported in Figure 2.
Also, the relative errors of the algorithms are shown in Figure 3. The
numerical results presented in Figures 2 and 3, show that Algorithm 7 is much
faster than Algorithm 6. We see that the proposed algorithms scale well to the
dimension of the data tensors. The accuracy achieved by the proposed
algorithms was also almost the same and even better than the baseline
algorithm (Algorithm 5), so this indicates the better performance and
efficiency of the proposed algorithms.

Figure 2: Example 1. Running time comparison of different algorithms.

Figure 3: Example 1. Relative error comparison of the algorithms

##  7 Conclusion and future works

In this paper, we proposed two fast randomized algorithms to compute the
generalized t-SVD (GTSVD) of tensors based on the tubal product (t-product).
Given two third-order tensors, the random projection technique is first used
to compute two small tensor sketches of the given tensors capturing the most
important ranges of them. Then from the small sketches, we recovered the
GTSVDs of the original data tensor from the GTSVD of the small tensor
sketches, which are easier to analyze. The computer simulations were conducted
convincing the feasibility and applicability of the proposed randomized
algorithm. The error analysis of the proposed algorithms using the power
iteration needs to be investigated and this will be our future research. We
plan to also develop randomized algorithms for the computation of the GTSVD to
be applicable for steaming data tensors, which arises in real-world
applications. The generalization of the proposed algorithm to higher order
tensors is our ongoing research work.

## References

  * [1] G. H. Golub, C. F. Van Loan, Matrix computations, JHU press, 2013.
  * [2] C. F. Van Loan, Generalizing the singular value decomposition, SIAM Journal on numerical Analysis 13 (1) (1976) 76–83.
  * [3] C. C. Paige, M. A. Saunders, Towards a generalized singular value decomposition, SIAM Journal on Numerical Analysis 18 (3) (1981) 398–405.
  * [4] P. C. Hansen, Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion, SIAM, 1998.
  * [5] O. Alter, P. O. Brown, D. Botstein, Generalized singular value decomposition for comparative analysis of genome-scale expression data sets of two different organisms, Proceedings of the National Academy of Sciences 100 (6) (2003) 3351–3356.
  * [6] S. P. Ponnapalli, M. A. Saunders, C. F. Van Loan, O. Alter, A higher-order generalized singular value decomposition for comparison of global mrna expression from multiple organisms, PloS one 6 (12) (2011) e28072.
  * [7] N. Halko, P.-G. Martinsson, J. A. Tropp, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM review 53 (2) (2011) 217–288.
  * [8] W. Wei, H. Zhang, X. Yang, X. Chen, Randomized generalized singular value decomposition, communications on Applied mathematics and computation 3 (1) (2021) 137–156.
  * [9] A. K. Saibaba, J. Hart, B. van Bloemen Waanders, Randomized algorithms for generalized singular value decomposition with application to sensitivity analysis, Numerical Linear Algebra with Applications 28 (4) (2021) e2364.
  * [10] S. Ahmadi-Asl, A. Cichocki, A. H. Phan, M. G. Asante-Mensah, M. M. Ghazani, T. Tanaka, I. Oseledets, Randomized algorithms for fast computation of low rank tensor ring model, Machine Learning: Science and Technology 2 (1) (2020) 011001\.
  * [11] S. Ahmadi-Asl, S. Abukhovich, M. G. Asante-Mensah, A. Cichocki, A. H. Phan, T. Tanaka, I. Oseledets, Randomized algorithms for computation of Tucker decomposition and higher order SVD (HOSVD), IEEE Access 9 (2021) 28684–28706.
  * [12] M. Che, Y. Wei, Randomized algorithms for the approximations of Tucker and the tensor train decompositions, Advances in Computational Mathematics 45 (1) (2019) 395–428.
  * [13] S. Ahmadi-Asl, A randomized algorithm for tensor singular value decomposition using an arbitrary number of passes, arXiv preprint arXiv:2207.12542 (2022).
  * [14] M. E. Kilmer, K. Braman, N. Hao, R. C. Hoover, Third-order tensors as operators on matrices: A theoretical and computational framework with applications in imaging, SIAM Journal on Matrix Analysis and Applications 34 (1) (2013) 148–172.
  * [15] M. E. Kilmer, C. D. Martin, Factorization strategies for third-order tensors, Linear Algebra and its Applications 435 (3) (2011) 641–658.
  * [16] L. R. Tucker, et al., The extension of factor analysis to three-dimensional matrices, Contributions to mathematical psychology 110119 (1964).
  * [17] L. R. Tucker, Some mathematical notes on three-mode factor analysis, Psychometrika 31 (3) (1966) 279–311.
  * [18] F. L. Hitchcock, The expression of a tensor or a polyadic as a sum of products, Journal of Mathematics and Physics 6 (1-4) (1927) 164–189.
  * [19] F. L. Hitchcock, Multiple invariants and generalized rank of a p-way matrix or tensor, Journal of Mathematics and Physics 7 (1-4) (1928) 39–79.
  * [20] Z.-H. He, M. K. Ng, C. Zeng, Generalized singular value decompositions for tensors and their applications., Numerical Mathematics: Theory, Methods & Applications 14 (3) (2021).
  * [21] Y. Zhang, X. Guo, P. Xie, Z. Cao, CS decomposition and GSVD for tensors based on the t-product, arXiv preprint arXiv:2106.16073 (2021).
  * [22] A. Cichocki, N. Lee, I. Oseledets, A.-H. Phan, Q. Zhao, D. P. Mandic, et al., Tensor networks for dimensionality reduction and large-scale optimization: Part 1 low-rank tensor decompositions, Foundations and Trends® in Machine Learning 9 (4-5) (2016) 249–429.
  * [23] C. Lu, J. Feng, Y. Chen, W. Liu, Z. Lin, S. Yan, Tensor robust principal component analysis with a new tensor nuclear norm, IEEE transactions on pattern analysis and machine intelligence 42 (4) (2019) 925–938.
  * [24] C. D. Martin, R. Shafer, B. LaRue, An order-p tensor factorization with applications in imaging, SIAM Journal on Scientific Computing 35 (1) (2013) A474–A490.
  * [25] S. Ahmadi-Asl, An efficient randomized fixed-precision algorithm for tensor singular value decomposition, Communications on Applied Mathematics and Computation (2022) 1–20.
  * [26] W. Yu, Y. Gu, Y. Li, Efficient randomized algorithms for the fixed-precision low-rank matrix approximation, SIAM Journal on Matrix Analysis and Applications 39 (3) (2018) 1339–1359.
  * [27] U. O. Ugwu, L. Reichel, Tensor regularization by truncated iteration: a comparison of some solution methods for large-scale linear discrete ill-posed problem with a t-product, arXiv preprint arXiv:2110.02485 (2021).
  * [28] U. O. Ugwu, Viterative tensor factorization based on Krylov subspace-type methods with applications to image processing, Kent State University, PhD Thesis,, 2021.