# Theoretical Modeling and Bio-inspired Trajectory Optimization of A Multiple-
locomotion Origami Robot

Keqi Zhu1,∗, Haotian Guo1,∗, Wei Yu1, Hassen Nigatu1,2, Tong Li3, Huixu
Dong1,2,† $*$ Keqi Zhu and Haotian Guo contributed equally to this work 1
Robot Perception and Grasp Laboratory(Grasp Lab), Zhejiang University,
Hangzhou 310058, China. Corresponding author Huixu Dong, e-mail:
<EMAIL_ADDRESS>2 Robotics Institute of Zhejiang University(robotics
research center of yuyao). 3 Department of Sports Science, Zhejiang
University, 310058, Hangzhou, ChinaThis work was supported by the startup fund
of "New Hundredd-Talent Program" of Zhejiang University under Grant
109100+194422R3/008 and Grant 109100$*$1942222R3/009 from Huixu Dong.

###### Abstract

Recent research on mobile robots has focused on increasing their adaptability
to unpredictable and unstructured environments using soft materials and
structures. However, the determination of key design parameters and control
over these compliant robots are predominantly iterated through experiments,
lacking a solid theoretical foundation. To improve their efficiency, this
paper aims to provide mathematics modeling over two locomotion, crawling and
swimming. Specifically, a dynamic model is first devised to reveal the
influence of the contact surfaces’ frictional coefficients on displacements in
different motion phases. Besides, a swimming kinematics model is provided
using coordinate transformation, based on which, we further develop an
algorithm that systematically plans human-like swimming gaits, with maximum
thrust obtained. The proposed algorithm is highly generalizable and has the
potential to be applied in other soft robots with multiple joints. Simulation
experiments have been conducted to illustrate the effectiveness of the
proposed modeling.

##  I Introduction

With the burgeoning advancements in materials science and structural science,
soft robotics has emerged as a field of significant interest and extensive
research, notably when contrasted with traditional robots, which are often
characterized by their intricate and sizable electromechanical systems[1].
These contemporary soft robots are distinguished by their apparent benefits,
namely their lightweight, compact construction, and the cost-effectiveness of
their mass fabrication capabilities. Nevertheless, the efficacy of these
robots within various environments is intrinsically linked to the optimization
of their structural design parameters and proper control strategies.

The design of the soft robot plays a decisive role in determining the upper
limit of the robot’s performance. Notably, Li et al. designed a soft robot
through both systematic experiments and theoretical analyses, which can be
actuated successfully in a field test in the Mariana Trench down to a depth of
10,900 meters[2]. However, many studies focus on the bio-inspiration analysis
and function realization of soft robots, and few of them further consider the
optimization of structural parameters to improve motion efficiency [3, 4].
Using the black box method, Morzadec et al. focused on shape optimization,
which can skip the complex modeling process but is time-consuming[5]. Vikas et
al. proposed a cost function to optimize the key parameters, without the
modeling[6]. Furthermore, a subset of research has concentrated on the
application of topological approaches[7] and differentiable simulation
techniques[8] for the optimization of design parameters. The modeling
processes employed by these methods are inherently more complex, and their
transferability to other robots appears to be somewhat limited. Currently,
despite the rapid development of soft robot design, there still exists a lack
of how to efficiently combine the structural characteristics to improve
robots’ performance.

Figure 1: The diagram of the origami bio-inspired robot, retrieved from [9].

The performance of the soft robots relies not only on their structural design
but heavily on their control strategy or planned gaits. Usually, this is
obtained via bio-inspiration analysis, due to the compliant and complex nature
of these soft robots. For example, Tang et al. developed a frog’s webbed feet-
inspired swimming robot, with the gait obtained experimentally [10]. This is
true of Yang’s research, where they determined the gaits through repetitive
trials [11]. These experiment-based design optimizations are time-consuming
and often fail to obtain the optimal gaits that make full use of the provided
robot design. Notable works include Li [12] and Kim [13]. In the former,
reinforcement learning is applied to enable a snake robot to master two
swimming patterns, either with the fastest swimming speed or with the least
power consumption. However, a simulation platform is needed to collect massive
data for the training. In the latter, a systematical dynamics model of a
multimodal underwater robot is established, enabling the robots to delicately
walk, swim, and actively control buoyancy. However, their methods are too
complex to be generalized to more soft robots. Providing an efficient method
that is compatible with a wider variety of soft robots is challenging.

Origami-based soft robots, due to the constraints introduced with the origami
structures, facilitate both the design and control of soft robots for various
applications and thus attract researchers’ attention [14, 15]. Without loss of
generality, the purpose of this study is to provide a mathematical analysis of
the physical model of such a robot in terms of crawling and swimming
locomotion, taking the multi-locomotion origami robot in [9] as an example (as
seen in Fig.1), to present design optimization and a generalizable gait
planning algorithm, facilitating more efficient crawling and swimming motions.
The proposed origami robot has octagon-origami structures as the body designed
for crawling and the twisted tower origami-based arm for swimming. Especially,
We qualitatively explore the effects of the proposed robot’s design parameters
on multiple-locomotion, which provides some insights for constructing a
multiple-locomotion origami robot, and an algorithm is proposed, providing
hints for how human-like optimized swimming gaits reasoning. We analytically
test the proposed robot with an excellent qualitative agreement with the
locomotion performance and the theoretical models. The contributions of the
paper are summarized below

  1. 1.

We construct theoretical modeling regarding the geometric configurations of
robots consisting of two origami patterns, octagon origami and twisted tower
origami, for crawling and swimming.

  2. 2.

We integrate the periodic dynamic model of the robots with the design
parameters and motion variables, enabling the quantification of the influence
of surface coefficient on crawling motion and providing hints for soft robot
design.

  3. 3.

An heuristic-based algorithm has been developed to optimize swimming gaits and
has the potential to be generalized to a variety of multi-joint articulated
structures consisting of soft actuation elements.

##  II Methodology

Taking the multi-locomotion origami robot in our previous work as an example
[9], in this section, we construct the mathematical models regarding the body
and arms to optimize their key design parameters and gaits.

###  II-A Modeling of Crawling Locomotion

####  II-A1 Force Analysis and Dynamic model

The body of the proposed robots has an octagon-origami structure(as shown in
Fig. 2-A) integrated with robot legs, which converts body deformation to
linear motion to achieve a crawling gait. It mimics an inchworm with
anisotropic friction, which has been adopted in other origami-inspired robots
with similar gaits [16]. The octagon-origami structure here keeps the same
functionality and control schemes as [17]. The bending deformation of octagon-
origami structures (robotic body) rises once the air pressure is exacted on
these two octagon-origami structures for producing bending moments, which
causes changes in interaction forces. Thus, the accurate actuation moment due
to the air pressure is equivalent to two actuation moments at the ends of
octagon-origami structures following the general beam theory, as seen in Fig.
2.

Here, the robot’s forward motion is modeled into two phases. Assuming a
generic time instant $t$, we define the rolling from the front leg as the
first half of a cycle (at time step $t-1$) and that from the rear leg the
second half (at time step $t$) [18]. Then, a discrete dynamic equation is
developed following Newton’s Second Law, namely

Figure 2: The origami sheet for the octagon-origami structure(B). The black
lines and the red dot lines represent mountain and valley folding patterns,
respectively in (A). The force sketches of crawling captured at the time
circle $t-1$ for the first half (B) and t for the second half(C). The sign ‘+’
represents the positive direction. $h$ denotes the height for the specified
leg and $d$ is the horizontal distance between the reference point and the
front foot for the given leg, respectively; $M_{a}$ is the equivalent
actuation torque. The middle points (red points)of the ends of the octagon-
origami structure are considered reference points.

$\begin{matrix}\begin{aligned} m\ddot{x}_{t-1}&=F_{F,t-1}^{f}-F_{F,t-1}^{r}\\\
m\ddot{y}_{t-1}&=F_{R,t-1}^{f}+F_{R,t-1}^{r}-m\cdot g\\\
m\ddot{x}_{t}&=F_{F,t}^{f}-F_{F,t}^{r}\\\
m\ddot{y}_{t}&=F_{R,t}^{f}+F_{R,t}^{r}-m\cdot g\end{aligned}\end{matrix}$ (1)

With $F_{F}$=$\mu F_{R}$, where $m$ denotes the mass of the robot; $\mu$ is
the friction coefficient; $g$ is the gravitational acceleration; $\ddot{x}$
and $\ddot{y}$ represent accelerations along the $x$-axis and $y$-axis,
respectively.

Combining the above equations, the ground reaction forces ($F_{R}$) during two
phases of the cycle can be expressed as

$\displaystyle F_{R,t-1}^{f}$
$\displaystyle=\frac{m}{2}\left(\frac{1}{\mu}\ddot{x}_{t-1}+\ddot{y}_{t-1}+g\right)$
(2) $\displaystyle F_{R,t-1}^{r}$
$\displaystyle=\frac{m}{2}\left(-\frac{1}{\mu}\ddot{x}_{t-1}+\ddot{y}_{t-1}+g\right)$
$\displaystyle F_{R,t}^{f}$
$\displaystyle=\frac{m}{2}\left(-\frac{2}{\mu}\ddot{x}_{t}+\ddot{y}_{t}+g\right)$
$\displaystyle F_{R,t}^{r}$
$\displaystyle=\frac{m}{2}\left(\frac{2}{\mu}\ddot{x}_{t}+\ddot{y}_{t}+g\right)$

Note, assuming a short time interval $\Delta T$, the accelerations are defined
as

$\ddot{l}_{n}=\frac{l_{n+1}-2\cdot l_{n}+l_{n-1}}{\bigtriangleup
t^{2}},l\in\left\\{x,y\right\\}$ (3)

Figure 3: The geometrical and crawling sketches of the proposed robot. (A) The
leg configuration with dimensions $h^{i},H^{i},b^{i},r^{i}$, $i=l$ for the
front leg and $i=r$ for the rear leg; (B), the projection lengths $m^{i}$
(along $x$-axis) and $n^{i}$ (along $y$-axis) directions; (C) The sketch of
the proposed robot with geometrical parameters; (D-E) The status simplified
sketches of the robot crawling in the first half and the second half.
$\alpha_{j}^{{i}^{\prime}}$ and $l_{j}^{i}$ symbolize the equivalent angle and
the leg’s length. $j$ being 0 or 1, denoting the first or second half. The
subscripts 0 and 1 of the circumference angle $\theta$, the chord length $a$,
and the chord inclination $\beta_{k}$ point to the sample time for a motion
cycle. The robotic height can be mainly determined by the equivalent angle
$\alpha_{j}^{{i}^{\prime}}$ and leg length $l_{j}^{f}$ of the front leg.

####  II-A2 Deformation Analysis of Crawling

As illustrated in Fig.3, the geometric deformation model is parameterized by
the design parameters of the legs, $h^{i},H^{i},b^{i},r^{i}$ ($i=l$ represents
the front leg; $i=r$ denotes the rear leg.) and the variables,
$\theta_{k},\alpha_{k}^{f},\alpha_{k}^{r}$ ($\theta_{k}$ is the circumference
angle of the bending body; $\alpha_{k}^{f}$ and $\alpha_{k}^{r}$ are the
inclinations of the front and rear legs at the time circle $k$, respectively).
We construct the model between design parameters and motion variables over the
cycles of the first and second halves. The equivalent length $l_{k}^{f}$ and
incline angle $\alpha_{k}^{{f}^{\prime}}$ of the front leg by replacing the
actual leg conFiguration for simplifying the calculation, with the projection
lengths $m^{f}$ and $n^{f}$ along the $x$-axis and $y$-axis, respectively.

$\begin{aligned}
m^{f}&=h^{f}\cos\alpha_{k}^{f}-\left(b^{f}-r^{f}\right)\sin\alpha_{k}^{f}+\left(H^{f}+r^{f}\right)\cos\alpha_{k}^{f}\\\
n^{f}&=h^{f}\sin\alpha_{k}^{f}+\left(b^{f}-r^{f}\right)\cos\alpha_{k}^{f}+\left(H^{f}+r^{f}\right)\sin\alpha_{k}^{f}\end{aligned}$

(4) $\displaystyle l_{k}^{f}$ $\displaystyle=\sqrt{m^{f^{2}}+n^{f^{2}}}$ (5)
$\displaystyle\alpha_{k}^{f^{\prime}}$
$\displaystyle=\arctan\left(\frac{n^{f}}{m^{f}}\right)$

where the subscript $k$ represents the time instant. Similarly, the
corresponding equivalent values of the rear leg are provided as

$\begin{aligned}
m^{r}&=h^{r}\cos\alpha_{k}^{r}-\left(b^{r}-r^{r}\right)\sin\alpha_{k}^{r}+\left(H^{r}+r^{r}\right)\cos\alpha_{k}^{r}\\\
n^{r}&=h^{r}\sin\alpha_{k}^{r}+\left(b^{r}-r^{r}\right)\cos\alpha_{k}^{r}+\left(H^{r}+r^{r}\right)\sin\alpha_{k}^{r}\end{aligned}$

(6) $\displaystyle l_{k}^{r}$ $\displaystyle=\sqrt{m^{r^{2}}+n^{r^{2}}}$ (7)
$\displaystyle\alpha_{k}^{r^{\prime}}$
$\displaystyle=\arctan\left(\frac{n^{r}}{m^{r}}\right)$

Fig.3 indicates that when combining Eq.(6) and Eq.(7), the mathematical model
of the chord length $a_{k}$, the chord inclination $\beta_{k}$ and the
circumference angle $\theta_{k}$ of the arc $L$ of the octagon-origami
structure can be built based on the geometrical constraints as follows:

$\beta_{k}=\arcsin\left(\frac{l_{k}^{f}\cos\alpha_{k}^{f^{\prime}}-l_{k}^{r}\cos\alpha_{k}^{r^{\prime}}}{\alpha_{k}}\right)$
(8) $a_{k}=\frac{L}{\theta_{k}}\sin\theta_{k}$ (9)
$\displaystyle\theta_{k}=\alpha_{k}^{f}+\beta_{k}$ (10)
$\displaystyle\theta_{k}=\alpha_{k}^{r}-\beta_{k}$

The equations above define the specific state of the robot. Combining them, we
have all other variables that can be expressed via $\theta_{k}$, which means
that the proposed robot has one degree of freedom for an octagon-origami
structure. $\Delta x_{(1)},\Delta y_{(1)}$ during the first half and $\Delta
x_{(2)},\Delta y_{(2)}$ for the second are the displacements of the center of
gravity along the $x$ and $y$ axes, respectively. The right superscript $i$ in
$\Delta x_{(1)}^{i}$ is $f$ as the front leg or $r$ as the rear leg, as shown
in Fig.3(D, E). Note that when the radius of the arc formed from the octagon-
origami structure is much larger than the chord of the arc, the arc length is
almost equivalent to the arc chord. Therefore, the displacements of the center
of the gravity are provided as

$\begin{matrix}\Delta
x_{(1)}+\frac{a_{1}}{2}\cdot\cos\beta_{1}=l_{1}^{f}\cdot\sin\alpha_{1}^{f^{\prime}}-l_{0}^{f}\cdot\sin\alpha_{0}^{f^{\prime}}+\frac{a_{0}}{2}\cdot\cos\beta_{0}+\Delta
x_{(1)}^{f}\\\ \Delta
x_{(2)}+\frac{a_{1}}{2}\cdot\cos\beta_{1}=l_{0}^{r}\cdot\sin\alpha_{0}^{r^{\prime}}-l_{1}^{r}\cdot\sin\alpha_{1}^{r^{\prime}}+\frac{a_{0}}{2}\cdot\cos\beta_{0}+\Delta
x_{(2)}^{r}\end{matrix}$

(11)

Due to rolling, the displacement of the leg is described as

$\begin{array}[]{l}\Delta
x_{(1)}^{f}=r^{f}\cdot\left(\alpha_{1}^{f^{\prime}}-\alpha_{0}^{f^{\prime}}\right)\\\
\Delta
x_{(2)}^{r}=r^{r}\cdot\left(\alpha_{0}^{r^{\prime}}-\alpha_{1}^{r^{\prime}}\right)\\\
\end{array}$ (12)

$\begin{matrix}\begin{aligned} &\Delta y_{(1)}=\Delta y_{(2)}=\\\
&\frac{l_{1}^{f}\cos\alpha_{1}^{f^{\prime}}-l_{1}^{r}\cos\alpha_{1}^{r^{\prime}}}{2}-\frac{l_{0}^{f}\cos\alpha_{0}^{f^{\prime}}-l_{0}^{r}\cos\alpha_{0}^{r^{\prime}}}{2}\end{aligned}\end{matrix}$

(13)

In terms of Eqs. (4-13), to arrive at a maximum displacement for each crawling
circle, the $\Delta x_{(1)}$ and $\Delta x_{(s)}$ should be enlarged.
Moreover, $\Delta y_{(1)}$ and $\Delta y_{(2)}$ tend to be minimized, which
enables the proposed robot to move stably.

###  II-B Modeling of Swimming Locomotion

####  II-B1 The Origami Twisted Tower

Inspired by human physiology, the proposed robot can realize human-like
swimming locomotion based on three serially connected origami twist towers
following the “generalized” Kresling pattern [19]. Fig.4-C presents the
geometry of the origami twisted tower, where $K,L,N,M$ stand for the vertexes
of the panel surface. The geometric parameters are defined: $a,b,l$ and $n$
represent the polygon side length, the panel side length, the diagonal valley-
crease length, and the number of polygon sides, respectively. In addition, the
angle $\delta$ is half the internal angle of the basal polygon, and $\lambda$,
the angle ratio, is a transformation metric between open and closed states.
With the rotation angle $\theta$ of the top polygon changing, the variation of
vertexes in circumradius $R$ results in scaling the height $H$ of the origami
structure, which can be expressed as:

$\vspace{-1.5mm}H=(b^{2}-2R^{2}(1-\cos\theta))^{\frac{1}{2}}$ (14)

The origami tower’s dimensions are critical for design efficiency and ease of
construction. Following Bhovad’s work [19], we manually set $a=b=15$ mm, with
a 95° angle, yielding an initial angle ($\theta_{int}\approx 7$°). To prevent
interference between the rotating origami towers and the robot’s body, we
ensured adequate space by setting the separation ($L_{s}$) slightly larger
than $a$, which is $16$ mm in this case, referring to Fig.4 and Fig.5 for
clarity. Additionally, to ensure the robotic arms’ effectiveness in pushing
water during an outreach movement, the first joint is positioned as low as the
design permits, maintaining a $10$ mm distance from the baseline. This careful
sizing ensures the robot’s arm movements are both functional and feasible.

Figure 4: The Three Stages of Breaststroke locomotion(A), The sketch of the
proposed robot(B); the front view(C-1), top-down view(C-2) and the
corresponding crease pattern(C-3) of the origami twisted tower, where the
solid black lines represent the mountain folds and the dashed black line
denote the valley folds, respectively; the coordinate transformations for the
first, second, third joint of the robotic arm(D); the top view of an origami
joint(E). 1, 2, and 3 represent the lifting, outreaching, and rotating
actions, and 4 indicates a flapping plate. $x_{i},y_{i}$ and $z_{i}$ denote
the coordinate reference frame of the $i$-th joint.

####  II-B2 Forward Kinematics of 3-DoF Twisted Towers

The forward kinematics of the 3-DoF origami twisted towers are established
following the DH conventions, to illustrate the motion of the robotic arm.
Considering that origami twisting tower output displacement and rotation
simultaneously, we divide each tower into a composite of revolute joint and
prismatic joint by adding a dummy coordinate. The red coordinate system in
Fig.4-E remains its orientation with the twisted tower’s deformation, serving
to effectively describe the tower’s axial elongation at the P joints.
Conversely, the blue coordinate system is attached to the tower’s surface to
articulate the rotational motion at the R joints. This dual-coordinate
approach allows for a precise representation of the tower’s complex motions.
Table I summarizes the D-H parameters. Thus, the transformation matrix from
joint $i$ to joint $i+1$ is:

$T_{i(i+1)}=\left[\begin{array}[]{cccc}\cos\theta^{i}&-\cos\alpha_{i}\sin\theta^{i}&\sin\alpha_{i}\sin\theta^{i}&a_{i}\cos\theta^{i}\\\
\sin\theta^{i}&\cos\alpha_{i}\cos\theta^{i}&-\sin\alpha_{i}\cos\theta^{i}&a_{i}\sin\theta^{i}\\\
0&\sin\alpha_{i}&\cos\alpha_{i}&d^{i}\\\ 0&0&0&1\end{array}\right]$

(15)

As shown in Fig.4-D, we have

$\centering\begin{matrix}\theta_{i}^{t}=\theta_{int}+\theta_{i}\\\
d_{0}=H_{0}+L_{s},d_{1}=H_{1}+L_{s},d_{2}=H_{2}+L_{p}\end{matrix}\@add@centering$
(16)

where subscript $i\in\left\\{1,2,3\right\\}$ denotes the $i$-th link of the
robotic arm; $\theta_{i}$ and $H_{i}$ represent the rotation angle and the
height with such a rotation angle; $L_{s}$ and $L_{p}$ denote the supporting
link length and the padding link length, as shown in Fig.4-E; $\theta_{i}^{t}$
and $d_{i}$ indicate the total rotation angle and the total length,
correspondingly. Substituting DH parameters into the Eq.15, the final $T$ can
be expressed with Eq.17:

(17)

where $c$ and $s$ represent the $cos$ and $sin$, respectively. Thus, the
reachable trajectories of the arm can be obtained.

Table I: The standard D-H parameters of the proposed robotic arm(Separation of
coupled variables)

$i$ | $\alpha_{i}$ | $a_{i}$ | $d_{i}$ | $\theta_{i}$ | offset
---|---|---|---|---|---
1 | 0 | 0 | $\widehat{H_{1}}$ | 0 | 0
${1}^{\prime}$ | $90^{\circ}$ | 0 | $L_{s}$ | $\widehat{\theta^{1}}$ | $180^{\circ}$
2 | 0 | $-90^{\circ}$ | $\widehat{H_{2}}/2$ | 0 | 0
${2}^{\prime}$ | $90^{\circ}$ | 0 | $L_{s}-a/2$ | $\widehat{\theta^{2}}$ | $90^{\circ}$
3 | 0 | 0 | $\widehat{H_{3}}$ | 0 | 0
${3}^{\prime}$ | 0 | 0 | $L_{p}$ | $\widehat{\theta^{3}}$ | $-80^{\circ}$

####  II-B3 Problem Description on Gait Optimization

The forward efficiency of the proposed robot relies not only on its structural
design but heavily on an optimized gait. Drag-based swimming, a common aquatic
locomotion mode, involves alternating protraction and retraction of paddle-
like limbs, with the power stroke exerting more force than the recovery stroke
[20]. Modulating the effective surface area of appendages has been proven an
effective solution. Ideally, Paddles should maximize their area during the
power stroke and minimize resistance during the non-propulsive phase by
reorienting into a stream-lined shape [21]. Taking inspiration from human
breaststroke swimming, an algorithm is proposed to provide hints for optimized
swimming gaits and can be generalized to various multi-joint articulated arms
consisting of soft actuation elements.

Technically in the breaststroke, arm stroke is divided into four phases,
namely Recovery, glide, out-sweep, and in-sweep, as seen in Fig. 6-A. These
phases involve constantly changing hand path width, depth, and pose to achieve
improved efficiency [22]. Since the Out-sweep does not contribute to
propulsion and thrust is primarily generated in the In-sweep [23], we simplify
the analysis to focus on the propulsive (In-sweep) and drag strokes (Out-sweep
and Recovery), aiming for maximum thrust and minimum drag, respectively. Here,
two key points are paid attention to, namely the start and end points of the
in-sweep phase. Specifically, the start is the point in the workspace that
spreads out the furthest, with the maximum $x$ value, and the end corresponds
to the point with the minimum $y$ value.

Referring to the blade element theory [24], the thrust force is defined as:

$F_{d}=\frac{1}{2}C_{d}\rho A_{p}V^{2}$ (18)

where $\rho$, $C_{d}$, $A_{p}$, and $V$ are fluid density, drag coefficient,
projected area, and relative velocity, respectively. To quantify the
efficiency along a given path, the forward strength is determined as the
integral of thrust force in its forward direction (y+), generated with a unit
surface area. At this point, the problem becomes a trajectory optimization
problem, where the optimization objective is to maximize and minimize the
volume of water pushed corresponding to the trajectory.

####  II-B4 Human-like Swimming Gaits Optimization

To solve this problem, we model the question as a graph theoretic problem in
the discredited space. Considering a uniform swimming motion, where $V$ is
constant, the thrust is linearly proportional to the projected area with the
appendages, as illustrated in Eq. 18. A rasterized joint-space point cloud is
obtained by isometric sampling of the joint angles, and then forward
kinematics is used to obtain the corresponding pose, including position and
orientation. The specific projected area is derived by projecting the normal
of the end plane with the y-axis. Thus, for each node, we calculate its
incremental thrust strength with its neighboring nodes in the joint space,
which is the multiplication of the difference in the y-axis with the averaged
projected areas. The goal is to find the path that maximizes and minimizes the
volume of thrust strength along the path. $A*$ algorithm is used for path
planning, which balances efficiency and optimality. Specifically, for the
thrust stroke, the heuristic function is defined as the largest projected area
multiplied by the difference in y-coordinate. This allows the heuristic
function to overestimate the thrust strength and ensure optimality so as to
provide a path with maximum thrust strength. Conversely, for the drag strokes,
the goal is to minimize the total dragging that impedes the motion. Thus, the
heuristic function is set as the multiplication of the smallest projected area
with the difference in y-coordinate values between the current point and the
target point. Namely,

$\begin{matrix}H_{thrust}(P_{n})=\mathop{max}\limits_{P_{n}\in
P}(A_{p}(Pn))|y_{n}-y_{goal}|\\\ H_{drag}(P_{n})=\mathop{min}\limits_{P_{n}\in
P}(A_{p}(Pn))|y_{n}-y_{goal}|\end{matrix}$ (19)

Figure 5: The dimensions of the proposed robot(A) and the first joint of the
proposed robot lifts arms in the water(B). The red dashed line indicates the
dimension baseline(A). The blue dash lines represent the water(B). A blown-up
view of the robotic foot is provided, which has two frictional surfaces.

##  III Results and Analysis

In this part, we analytically evaluate the proposed robot’s multiple-
locomotion capabilities based on the proposed modeling, including crawling
analysis and swimming trajectory optimization.

###  III-A Crawling Locomotion Friction Optimization

The built dynamic force model is applied to determine the suitable design for
achieving the optimal movement speed and stability of the robot within the
constraints of geometrical ranges (see Fig.5), with a maximum height of 85mm
and length of 110mm.

Figure 6: (A) The workspace of the 3-DOF origami-based arm and the obtained
optimal swimming arm gaits, including thrust stroke (in-sweep) and drag stroke
(out-sweep and recovery); (B) The joints trajectory of the obtained gaits. (C)
The relationship of nodal projected area (efficiency) along the y-axis

####  III-A1 The robotic foot design based on the simulation

The dynamic model is used for determining the surface materials to determine
the proposed legs’ surface materials by evaluating the ground reaction forces
during one cycle. The corresponding parameters have been set, namely
$g=9.8(m/s^{2})$, and $m=0.1(kg)$. For simplicity, we consider five discrete
values of friction coefficients ($\mu\in\left\\{0.1k\mid k=1,...,5\right\\}$)
to investigate their influence on crawling efficiency. Fig.7(A-1, A-2) and
(A-3, A-4) depict the magnitudes of the ground reaction forces experienced by
the front and rear feet at sequential time circles $t-1$ and $t$,
respectively, derived by Eqs. (2-13). Ground reaction forces are computed by
multiplying directly with the friction coefficient to ascertain the contact
frictional force between the foot and the ground, shown in Fig. 2. As
evidenced in Fig. 2 and Fig. 5-B, during the time circle $t-1$, the toes of
the rear foot and the pelma surface of the front foot are concurrently in
contact with the ground. To optimize forward locomotion at this juncture,
maximum the distance traveled in a time cycle(Eqs. (12,13)), it is crucial to
immobilize the front foot while propelling the rear foot forward by leveraging
the frictional interaction with the ground. The converse is applicable at time
circle $t$.

Consequently, a superior frictional coefficient is attributed to the plantar
surfaces of both the front and rear feet compared to the toe sections. This
differential in frictional attributes facilitates maximal forward progression
during the time circle $t-1$ and minimizes resistive effects at time circle
$t$.

Figure 7: The ground reaction forces $F_{R,t-1}^{f}$, $F_{R,t-1}^{r}$ for the
first half cycle(A-1, A-2) and $F_{R,t}^{f}$, $F_{R,t}^{r}$ for the second
half cycle (A-3, A-4) in terms of the circle angles $\theta_{t-1}$,
$\theta_{t}$ of arcs actuated by the control system at the time circles $t-1$
and $t$, and (B) shows the friction, respectively.

###  III-B Swimming Locomotion Trajectory Optimization

####  III-B1 The optimization results and simulation

Fig.6-A presents the workspace of the 3-DOF origami-based arm and obtained
optimal swimming arm gaits, including thrust stroke (in-sweep) and drag stroke
(out-sweep and recovery), with the color corresponding to the size of the unit
projected area at the given point. It can be seen that the proposed two-state
planned trajectory aligns well with the arm trajectory of human breaststroke
swimming (Fig. 4-A).

Due to the mechanical characteristics and the spatial alignment of the three
origami twisted towers, the resultant workspace is a piece of a thick
spherical shell. It should be noted that at the beginning of the recovery
stroke, the planned trajectory passes through the workspace point cloud, with
the rapid decrease of $\theta_{3}$ as seen in Fig. 6-B. This can be ascribed
to the large difference in the projected area at the interior and exterior of
the workspace in Cartesian space. The conduct of the trajectory search in the
joint space rather than Cartesian space provides trajectories that exhibit
superior continuity and smoothness, avoiding over-quick deformation that is
not feasible for fragile soft structures. With systematical gait planning, the
proposed robot can rapidly adjust the orientation of the end at the junction
of the two stages, shifting the projected area from a maximum case to the
minimum state in the shortest distance. Thus we alleviate the influence of
dragging in the swimming cycle, demonstrating the effectiveness of our
proposed algorithm. Fig.6-C presents the relationship of nodal projected area
(efficiency) along the y-axis in the workspace. The area enclosed by a
specific curve is equivalent to the thrust/drag produced by that gait.

To better understand the posture of the robot while swimming, we built a
simulation model of the robot in Matlab. Each side of the twisted tower was
accurately modeled to obtain a more realistic swimming state. Using the three-
degree-of-freedom robotic arm joint rotation angle we obtained through the A*
trajectory search algorithm, we simulated the posture of the robot’s
breaststroke throughout the cycle. From Fig.8-A to Fig.8-C is the entire
insweep stage, which is responsible for efficiently giving the robot forward
force; Fig.8-C to Fig.8-D, Fig.8-D to Fig.8-B are the recovery stage
respectively. and outweep stages, our algorithm ensures that in these two
stages, the force hindering the robot’s progress is minimal. What is worth
noting is the change in posture from Fig.8-B to Fig.8-A. From the end of the
previous cycle to the beginning of the next cycle, the third joint will rotate
independently, converting the end of the arm from a low-resistance water-
cutting posture to a high-resistance The water-pushing posture is more in line
with bionics. Fig.8-E, F show the trajectory left by the robot during its
advancement under multiple action cycles. Compared with our previous
research[9], the swimming efficiency has been greatly improved, and it is more
like a breaststroke posture.

###  III-C Discussion

The physical performance of the multiple-locomotion capabilities based on the
proposed modeling has been verified and presented in our previous work [9], as
seen in Fig. 9. With the optimized parameters and gaits, the proposed robot is
able to perform stable crawling and swimming locomotion.

We have provided the corresponding comparisons with other works. The parameter
optimization method we proposed is based on simplified beam theory with linear
deformation, while Soft materials such as rubber, silicone, etc. often have
complex nonlinear relationships. Compared with the topology-based methods with
finite elements analysis (FEA) illustrated in [7], we only obtain the
qualitative relationship among the design parameters. Besides, in contrast to
the dynamic modeling proposed in [13], the proposed robot has shortcomings in
the capability of underwater locomotion and velocity control for buoyancy and
turning. The current method plans the gaits in the position map by assuming a
constant forward velocity, neglecting the positive influence of velocity on
the dragging force produced by water. However, this simple strategy can be
easily generalized to a wider range of soft robots to improve their
efficiency. Meanwhile, although the body is designed for crawling, swimming
efficiency can be further improved through arm-leg coordination by considering
the contribution and impedance of the legs.

Figure 8: The snapshots of the robot swimming locomotion, (A)end of out-sweep;
(B) beginning of in-sweep; (C) beginning of recover; (D) beginning of out-
sweep; (E) and (F) demonstrate the trajectory of the end of the arm as the
robot swims forward.

Figure 9: The snapshots of the robot experiment in physical scenarios[9] (A)
demonstrate swimming locomotion, and (B) for crawling locomotion

##  IV CONCLUSIONS

In this work, we constructed mathematical models to provide the theoretical
foundation for the design optimization and control of soft origami robots on
crawling and swimming locomotion. By combining the periodic dynamic model of
the robots with the design parameters and motion variables, we qualitatively
reveal the influence of the surface frictional coefficient on crawling motion,
and provide hints for crawling-based robot design. Besides, a swimming
kinematics model, together with a heuristic-based algorithm, has been
developed to optimize the swimming trajectory with the least effort. Both
simulations and experiments were carried out to illustrate the effectiveness
of the proposed modeling and gaits control strategies, which allow the
proposed robot to successfully imitate an inchworm that crawls on the ground
and humans that swim in the water, respectively. Future research will be
focused on improving the capability of underwater locomotion for such robots.
We will continue exploring the dynamic modeling of the proposed robot to
further optimize the design and research on arm-leg coordination for improved
performance.

## References

  * [1] Junzhi Yu et al.  “On a Bio-inspired Amphibious Robot Capable of Multimodal Motion”  In _IEEE/ASME Transactions on Mechatronics_ 17.5, 2012, pp. 847–856  DOI: 10.1109/TMECH.2011.2132732
  * [2] Guorui Li et al.  “Self-powered soft robot in the Mariana Trench”  In _Nature_ 591.7848 Nature Publishing Group UK London, 2021, pp. 66–71
  * [3] Qiuxuan Wu et al.  “A novel underwater bipedal walking soft robot bio-inspired by the coconut octopus”  In _Bioinspiration & Biomimetics_ 16.4 IOP Publishing, 2021, pp. 046007
  * [4] Zhiji Han, Zhijie Liu, Wei He and Guang Li  “Distributed parameter modeling and boundary control of an octopus tentacle-inspired soft robot”  In _IEEE Transactions on Control Systems Technology_ 30.3 IEEE, 2021, pp. 1244–1256
  * [5] Thomas Morzadec, Damien Marcha and Christian Duriez  “Toward shape optimization of soft robots”  In _2019 2nd IEEE International Conference on Soft Robotics (RoboSoft)_ , 2019, pp. 521–526  IEEE
  * [6] Vishesh Vikas et al.  “Design and locomotion control of a soft robot using friction manipulation and motor–tendon actuation”  In _IEEE Transactions on Robotics_ 32.4 IEEE, 2016, pp. 949–959
  * [7] Brandon Caasenbrood, Alexander Pogromsky and Henk Nijmeijer  “A computational design framework for pressure-driven soft robots through nonlinear topology optimization”  In _2020 3rd IEEE international conference on soft robotics (RoboSoft)_ , 2020, pp. 633–638  IEEE
  * [8] Moritz Bächer, Espen Knoop and Christian Schumacher  “Design and control of soft robots using differentiable simulation”  In _Current Robotics Reports_ 2.2 Springer, 2021, pp. 211–221
  * [9] Huixu Dong et al.  “Bioinspired Amphibious Origami Robot with Body Sensing for Multimodal Locomotion” PMID: 35671518  In _Soft Robotics_ 9.6, 2022, pp. 1198–1209  DOI: 10.1089/soro.2021.0118
  * [10] Yucheng Tang et al.  “A frog-inspired swimming robot based on dielectric elastomer actuators”  In _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2017, pp. 2403–2408  DOI: 10.1109/IROS.2017.8206054
  * [11] Y. Yang et al.  “3D-Printed Origami Actuators for a Multianimal-Inspired Soft Robot with Amphibious Locomotion and Tongue Hunting” Advance online publication  In _Soft robotics_ , 2024  DOI: 10.1089/soro.2023.0079
  * [12] Guanda Li, Jun Shintake and Mitsuhiro Hayashibe  “Soft-body dynamics induces energy efficiency in undulatory swimming: A deep learning study”  In _Frontiers in Robotics and AI_ 10 Frontiers Media SA, 2023, pp. 1102854
  * [13] Taesik Kim, Juhwan Kim and Son-Cheol Yu  “Development of Bioinspired Multimodal Underwater Robot “HERO-BLUE” for Walking, Swimming, and Crawling”  In _IEEE Transactions on Robotics_ 40, 2024, pp. 1421–1438  DOI: 10.1109/TRO.2024.3353040
  * [14] Sébastien D. Rivaz et al.  “Inverted and vertical climbing of a quadrupedal microrobot using electroadhesion”  In _Science Robotics_ 3.25, 2018, pp. eaau3038  DOI: 10.1126/scirobotics.aau3038
  * [15] Daniela Rus and Michael T. Tolley  “Design, fabrication and control of origami robots”  In _Nature Reviews Materials_ 3.6, 2018, pp. 101–112  DOI: 10.1038/s41578-018-0009-8
  * [16] Je-Sung Koh and Kyu-Jin Cho  “Omega-Shaped Inchworm-Inspired Crawling Robot With Large-Index-and-Pitch (LIP) SMA Spring Actuators”  In _IEEE/ASME Transactions on Mechatronics_ 18.2, 2013, pp. 419–429  DOI: 10.1109/TMECH.2012.2211033
  * [17] Laura Paez, Maritza Granados and Kamilo Melo  “Conceptual design of a modular snake origami robot”  In _2013 IEEE International Symposium on Safety, Security, and Rescue Robotics (SSRR)_ , 2013, pp. 1–2  DOI: 10.1109/SSRR.2013.6719379
  * [18] N. Lobontiu, M. Goldfarb and E. Garcia  “A piezoelectric-driven inchworm locomotion device”  In _Mechanism and Machine Theory_ 36.4, 2001, pp. 425–443  DOI: https://doi.org/10.1016/S0094-114X(00)00056-2
  * [19] Priyanka Bhovad and Suyi Li  “Using Multi-Stable Origami Mechanism for Peristaltic Gait Generation: A Case Study”  In _Volume 5B: 42nd Mechanisms and Robotics Conference_ Quebec City, Quebec, Canada: American Society of Mechanical Engineers, 2018, pp. V05BT07A061  DOI: 10.1115/DETC2018-85932
  * [20] Bokeon Kwak, Dongyoung Lee and Joonbum Bae  “Efficient Drag-based Swimming using Articulated Legs with Micro Hair Arrays Inspired by a Water Beetle”  In _2019 2nd IEEE International Conference on Soft Robotics (RoboSoft)_ , 2019, pp. 795–800  DOI: 10.1109/ROBOSOFT.2019.8722728
  * [21] Naohiko Nishijima Hideki Takagi and Barry Wilson  “Swimming” PMID: 15079985  In _Sports Biomechanics_ 3.1 Routledge, 2004, pp. 15–27  DOI: 10.1080/14763140408522827
  * [22] Vassilios Gourgoulis and Thomas Nikodelis  “Comparison of the arm-stroke kinematics between maximal and sub-maximal breaststroke swimming using discrete data and time series analysis”  In _Journal of Biomechanics_ 142, 2022, pp. 111255  DOI: https://doi.org/10.1016/j.jbiomech.2022.111255
  * [23] E.W. Maglischo  “Swimming Even Faster”  Mayfield Publishing Company, 1993  URL: https://books.google.com/books?id=v4BZAAAAYAAJ
  * [24] Julianna M Gal and RW Blake  “Biomechanics of frog swimming: II. Mechanics of the limb-beat cycle in Hymenochirus Boettgeri”  In _Journal of experimental biology_ 138.1 The Company of Biologists Ltd, 1988, pp. 413–429


# Co-Clustering for Multitask Learning

Keerthiram Murugesan  Jaime Carbonell  Yiming Yang

###### Abstract

This paper presents a new multitask learning framework that learns a shared
representation among the tasks, incorporating both task and feature clusters.
The jointly-induced clusters yield a shared latent subspace where task
relationships are learned more effectively and more generally than in state-
of-the-art multitask learning methods. The proposed general framework enables
the derivation of more specific or restricted state-of-the-art multitask
methods. The paper also proposes a highly-scalable multitask learning
algorithm, based on the new framework, using conjugate gradient descent and
generalized Sylvester equations. Experimental results on synthetic and
benchmark datasets show that the proposed method systematically outperforms
several state-of-the-art multitask learning methods.

Multitask learning, Co-clustering, Matrix Factorization

##  1 Introduction

Multitask learning leverages shared structures among the tasks to jointly
build a better model for each task. Most existing work in multitask learning
focuses on how to take advantage of task similarities, either by learning the
relationship between the tasks via cross-task regularization techniques (Zhang
& Yeung, 2014; Zhang & Schneider, 2010; Rothman et al., 2010; Xue et al.,
2007) or by learning a shared feature representation across all the tasks,
leveraging low-dimensional subspaces in the feature space (Argyriou et al.,
2008; Jalali et al., 2010; Liu et al., 2009; Swirszcz & Lozano, 2012).
Learning task relationships has been shown beneficial in (positive and
negative) transfer of knowledge from information-rich tasks to information-
poor tasks (Zhang & Yeung, 2014), whereas the shared feature representation
has been shown to perform well when each task has a limited number of training
instances (observations) compared to the total number across all tasks
(Argyriou et al., 2008). Existing research in multitask learning considers
either the first approach and learns a task relationship matrix in addition to
the task parameters, or relies on the latter approach and learns a shared
latent feature representation from the task parameters. To the best of our
knowledge, there is no prior work that utilizes both principles jointly for
multitask learning. In this paper, we propose a new approach that learns a
shared feature representation along with the task relationship matrix jointly
to combine the advantages of both principles into a general multitask learning
framework.

Early work on latent shared representation includes (Zhang et al., 2005),
which proposes a model based on Independent Component Analysis (ICA) for
learning multiple related tasks. The task parameters are assumed to be
generated from independent sources. (Argyriou et al., 2008) consider sparse
representations common across many learning tasks. Similar in spirit to PCA
for unsupervised tasks, their approach learns a low dimensional representation
of the observations (Ding & He, 2004). More recently, (Kumar & Daume, 2012)
assume that relationships among tasks are sparse to enforce that each observed
task is obtained from only a few of the latent features, and from there learn
the overlapping group structure among the tasks. (Crammer & Mansour, 2012)
propose a K-means-like procedure that simultaneously clustering different
tasks and learning a small pool of $m\ll T$ shared models. Specifically, each
task is free to choose a model from the pool that better classifies its own
data, and each model is learned from pooling together all the training data
that belongs to the same cluster. (Barzilai & Crammer, 2015) propose a similar
approach that clusters the $T$ tasks into $K$ task-clusters with hard
assignments.

These methods compute the factorization of the task weight matrix to learn the
shared feature representation and the task structure. This matrix
factorization induces the simultaneous clustering of both the tasks and the
features in the $K$-dimensional latent subspace (Li & Ding, 2006). One of the
major disadvantages of this assumption is that it restricts the model to
define both the tasks and the features to have same number of clusters. For
example, in the case of sentiment analysis, where each task belongs to a
certain domain or a product category such as books, automobiles, etc., and
each feature is simply a word from the vocabulary of the product reviews.
Clearly, assuming both the features and the tasks have same number of clusters
is an unjustified assumption, as the number of feature clusters are typically
more than the number of task clusters, but the latter increase more than the
former, as new products are introduced. Such a restrictive assumption may (and
often does) hurt the performance of the model.

Unlike in the previous work, our proposed approach provides a flexible way to
cluster both the tasks and the features. We introduce an additional degree of
freedom that allows the number of task clusters to differ from the number of
features clusters (Ding et al., 2006; Wang et al., 2011). In addition, our
proposed models learns both the task relationship matrix and the feature
relationship matrix along with the co-clustering of both the tasks and the
features (Gu & Zhou, 2009; Sindhwani et al., 2009). Our proposed approach is
closely related to Output Kernel Learning (OKL) where we learn the kernel
between the components of the output vector for problems such as multi-output
learning, multitask learning, etc (Dinuzzo et al., 2011; Sindhwani et al.,
2013). The key disadvantage of OKL is that it requires the computation of
kernel matrix between every pair of instances from all the tasks. This results
in scalability constraint especially when the number of tasks/features is
large (Weinberger et al., 2009). Our proposed models achieve the similar
effect by learning a shared feature representation common across the tasks.

A key challenge in factoring with the extra degree of freedom is optimizing
the resulting objective function. Previous work on co-clustering for multitask
learning requires strong assumptions on the task parameters. (Zhong & Kwok,
2012) or not scalable to large-scale applications (Xu et al., 2015). We
propose an efficient algorithm that scales well to large-scale multitask
learning and utilizes the structure of the objective function to learn the
factorized task parameters. We formulate the learning of latent variables in
terms of a generalized Sylvester equation which can be efficiently solved
using the conjugate gradient descent algorithm. We start from the mathematical
background and then motivate our approach in Section 2. Then we introduce our
proposed models and their learning procedures in Section 3. Section 4 reports
the empirical analysis of our proposed models and shows that learning both the
task clusters and the feature clusters along with the task parameters gives
significant improvements compared to the state-of-the-art baselines in
multitask learning.

##  2 Background

Suppose we have $T$ tasks and
$\mathcal{D}_{t}=\\{\mathbf{X}_{t},Y_{t}\\}=\\{(x_{ti},y_{ti}):i=1,2,...,N_{t}\\}$
is the training set for each task $t=\\{1,2,\ldots,T\\}$. Let $W_{t}$
represent the weight vector for a task indexed by $t$. These task weight
vectors are stacked as columns of a matrix $\mathbf{W}$, which is of size
$P\times T$ , with $P$ being the feature dimension. Traditional multitask
learning imposes additional assumptions on $\mathbf{W}$ such as low-rank,
$\ell_{1}$ norm, $\ell_{2,1}$ norm, etc to leverage the shared characteristics
among the tasks. In this paper, we consider a similar assumption based on the
factorization of the task weight matrix $\mathbf{W}$.

In factored models, we decompose the weight matrix $\mathbf{W}$ as
$\mathbf{\mathbf{missing}}F\mathbf{G}^{\top}$, where $\mathbf{F}$ can be
interpreted as a feature cluster matrix of size $P\times K$ with $K$ feature
clusters and, similarly, $\mathbf{G}$ as a task cluster matrix of size
$T\times K$ with $K$ task clusters. If we consider squared error losses for
all the tasks, then the objective function for learning $\mathbf{F}$ and
$\mathbf{G}$ can be given as follows:

$\displaystyle\operatorname*{arg\,min}_{\tiny\begin{subarray}{c}\mathbf{F}\in\mathcal{R}^{P\times
K}\\\ \mathbf{G}\in\mathcal{R}^{T\times K}\\\
\mathbf{F}\in\Gamma_{\mathcal{F}},\mathbf{G}\in\Gamma_{\mathcal{G}}\end{subarray}}$
$\displaystyle\sum_{t\in[T]}\lVert{Y_{t}-\mathbf{X}_{t}\mathbf{F}G_{t}^{\top}}\rVert_{2}^{2}+\mathcal{P}_{\lambda_{1}}(\mathbf{F})+\mathcal{P}_{\lambda_{2}}(\mathbf{G})$
(1)

In the above objective function, the latent feature representation is captured
by the matrix $\mathbf{F}$ and the grouping structure on the tasks is
determined by the matrix $\mathbf{G}$. The predictor $W_{t}$ for task $t$ can
then be computed from $\mathbf{F}G_{t}^{\top}$, where $G_{t}$ is $t^{th}$ row
of matrix $\mathbf{G}$. In the above objective function,
$\mathcal{P}_{\lambda_{1}}(\mathbf{F})$ is a regularization term that
penalizes the unknown matrix $\mathbf{F}$ with regularization parameter
$\lambda_{1}$. Similarly, $\mathcal{P}_{\lambda_{2}}(\mathbf{G})$ is a
regularization term that penalizes the unknown matrix $\mathbf{G}$ with
regularization parameter $\lambda_{2}$. $\Gamma_{\mathcal{F}}$ and
$\Gamma_{\mathcal{G}}$ are their corresponding constraint spaces. Without
these additional constraints on $\mathbf{F}$ and $\mathbf{G}$, the objective
function reduces to solving each task independently, since any task weight
matrix from $\mathbf{F}$ and $\mathbf{G}$ can also be attained by
$\mathbf{W}$.

Several assumptions can be enforced on these unknown factors $\mathbf{F}$ and
$\mathbf{G}$. Below we discuss some of the previous models that make some
well-known assumptions on $\mathbf{F}$ and $\mathbf{G}$ and can be written in
terms of the above objective function.

(1) Factored Multitask Learning (FMTL) (Amit et al., 2007) considers a squared
frobenius norm on both $\mathbf{F}$ and $\mathbf{G}$.

$\displaystyle\operatorname*{arg\,min}_{\tiny\begin{subarray}{c}\mathbf{F}\in\mathcal{R}^{P\times
K}\\\ \mathbf{G}\in\mathcal{R}^{T\times K}\end{subarray}}$
$\displaystyle\sum_{t\in[T]}\lVert{Y_{t}-\mathbf{X}_{t}\mathbf{F}G_{t}^{\top}}\rVert_{2}^{2}+\lambda_{1}\lVert{\mathbf{F}}\rVert_{F}^{2}+\lambda_{2}\lVert{\mathbf{G}}\rVert_{F}^{2}$
(2)

It can be shown that the above problem can equivalently written as the
multitask learning with trace norm constraint on the task weight matrix
$\mathbf{W}$.

(2) Multitask Feature Learning (MTFL) (Argyriou et al., 2008) assumes that the
matrix $\mathbf{G}$ learns sparse representations common across many tasks.
Similar in spirit to PCA for unsupervised tasks, MTFL learns a low dimensional
representation of the observations $\mathbf{X}_{t}$ for each task, using
$\mathbf{F}$ such that $\mathbf{F}\mathbf{F}^{\top}=\mathbf{I}_{p}$.

$\displaystyle\operatorname*{arg\,min}_{\begin{subarray}{c}\mathbf{F}\in\mathcal{R}^{P\times
K},\mathbf{G}\in\mathcal{R}^{T\times K}\\\
\mathbf{F}\mathbf{F}^{\top}=\mathbf{I}_{p}\end{subarray}}$
$\displaystyle\sum_{t\in[T]}\lVert{Y_{t}-\mathbf{X}_{t}\mathbf{F}G_{t}^{\top}}\rVert_{2}^{2}+\lambda\lVert{\mathbf{G}}\rVert_{2,1}^{2}$
(3)

where $K$ is usually set to $P$. It considers an $\ell_{2,1}$ norm on
$\mathbf{G}$ to force all the tasks to have a similar sparsity pattern such
that the tasks select the same latent features (columns of $\mathbf{F}$). It
is worth noting that the Equation 3 can be equivalently written as follows:

$\displaystyle\operatorname*{arg\,min}_{\begin{subarray}{c}\mathbf{W}\in\mathcal{R}^{P\times
T},\\\ \mathbf{\Sigma}\succeq 0\end{subarray}}$
$\displaystyle\sum_{t\in[T]}\lVert{Y_{t}-\mathbf{X}_{t}W_{t}}\rVert_{2}^{2}+\lambda
tr({\mathbf{W}^{\top}\mathbf{\Sigma}^{-1}\mathbf{W}})$ (4)

which then can be rewritten as multitask learning with a trace norm constraint
on the task weight matrix $\mathbf{W}$ as before.

(3) Group Overlap MTL (GO-MTL) (Kumar & Daume, 2012) assumes that the matrix
$\mathbf{G}$ is sparse to enforce that each observed task is obtained from
only a few of the latent features, indexed by the non-zero pattern of the
corresponding rows of the matrix $\mathbf{G}$.

$\displaystyle\operatorname*{arg\,min}_{\tiny\begin{subarray}{c}\mathbf{F}\in\mathcal{R}^{P\times
K}\\\ \mathbf{G}\in\mathcal{R}^{T\times K}\end{subarray}}$
$\displaystyle\sum_{t\in[T]}\lVert{Y_{t}-\mathbf{X}_{t}\mathbf{F}G_{t}^{\top}}\rVert_{2}^{2}+\lambda_{1}\lVert{\mathbf{F}}\rVert_{F}^{2}+\lambda_{2}\lVert{\mathbf{G}}\rVert_{1}^{1}$
(5)

The above objective function can be compared to dictionary learning where each
column of $\mathbf{F}$ is considered as a dictionary atom and each row of
$\mathbf{G}$ is considered as their corresponding sparse codes (Maurer et al.,
2013).

(4) Multitask Learning by Clustering (CMTL) (Barzilai & Crammer, 2015) assumes
that the $T$ tasks can be clustered into $K$ task-clusters with hard
assignment. For example, if the $k$th element of $G_{t}$ is one, and all other
elements of $G_{t}$ are zero, we say that task $t$ is associated with cluster
$k$.

$\displaystyle\operatorname*{arg\,min}_{\tiny\begin{subarray}{c}\mathbf{F}\in\mathcal{R}^{P\times
K},\mathbf{G}\in\mathcal{R}^{T\times K}\\\ G_{t}\in\\{0,1\\}^{K}\\\
\lVert{G_{t}}\rVert_{2}=1,\forall t\in[T]\end{subarray}}$
$\displaystyle\sum_{t\in[T]}\lVert{Y_{t}-\mathbf{X}_{t}\mathbf{F}G_{t}^{\top}}\rVert_{2}^{2}+\lambda_{1}\lVert{\mathbf{F}}\rVert_{F}^{2}$
(6)

The constraints $G_{t}\in\\{0,1\\}^{K},\lVert{G_{t}}\rVert_{2}=1$ ensure that
$\mathbf{G}$ is a proper clustering matrix. Since the above problem is
computationally expensive as it involves solving a combinatorial problem, the
constraint on $\mathbf{G}$ is relaxed as $G_{t}\in[0,1]^{K}$.

These four methods require the number of task clusters to be same as the
number of features clusters, which as mentioned earlier, is a restrictive
assumption that may and often does hurt performance. In addition, these
methods do not leverage the inherent relationship between the features (via
$\mathbf{F}$) and the relationship between the tasks (via $\mathbf{G}$). Note
that these objective functions are bi-convex problems where the optimization
is convex in $\mathbf{F}$ when fixing $\mathbf{G}$ and vice versa. We cannot
achieve globally optimal solution but one can show that algorithm reaches the
locally optimal solution in a fixed number of iterations.

##  3 Proposed Approach

###  3.1 BiFactor MTL

Existing models do not take into consideration both the relationship between
the tasks and the relationship between the features. Here we consider a more
general formulation that in addition to estimating the parameters $\mathbf{F}$
and $\mathbf{G}$, we learn their task relationship matrix $\mathbf{\Omega}$
and the feature relationship matrix $\mathbf{\Sigma}$. We call this framework
BiFactor multitask learning, following the factorization of the task
parameters $\mathbf{W}$ into two low-rank matrices $\mathbf{F}$ and
$\mathbf{G}$.

$\displaystyle\operatorname*{arg\,min}_{\tiny\begin{subarray}{c}\mathbf{F}\in\mathcal{R}^{P\times
K},\mathbf{G}\in\mathcal{R}^{T\times K}\\\ \mathbf{\Sigma}\succeq
0,\mathbf{\Omega}\succeq 0\end{subarray}}$
$\displaystyle\sum_{t\in[T]}\lVert{Y_{t}-\mathbf{X}_{t}\mathbf{F}G_{t}}\rVert_{2}^{2}$
(7)
$\displaystyle+\lambda_{1}tr({\mathbf{F}^{\top}\mathbf{\Sigma}^{-1}\mathbf{F}})+\lambda_{2}tr({\mathbf{G}^{\top}\mathbf{\Omega}^{-1}\mathbf{G}})$

In the above objective function, we consider
$\mathcal{P}_{\lambda_{1}}(\mathbf{F})$ and
$\mathcal{P}_{\lambda_{2}}(\mathbf{G})$ to learn task relationship and feature
relationship matrices $\mathbf{\Sigma}$ and $\mathbf{\Omega}$. The motivation
for these regularization terms is based on (Argyriou et al., 2008; Zhang &
Yeung, 2014) where they considered separately either the task relationship
matrix $\mathbf{\Omega}$ or the feature relationship matrix $\mathbf{\Sigma}$.
Note that the value of $K$ is typically set to value less than $\min(P,T)$.

It is easy to see that by setting the value of $\mathbf{G}$ to $I_{T}$
111identity matrix of size $T\times T$ (assuming that the rank $K$ is set to
$T$) , our objective function reduces to multitask feature learning ($MTFL$)
discussed in the previous section. Similarly, by setting the value of
$\mathbf{F}$ to $I_{P}$ 222identity matrix of size $P\times P$ (assuming that
the rank $K$ is set to $P$) , our objective function reduces to multitask
relationship learning ($MTRL$) (Zhang & Yeung, 2014). If we set
$\mathbf{\Omega}=I_{T}$ and $\mathbf{\Sigma}=I_{P}$, we obtain the factored
multitask learning setting ($FMTL$) defined in Equation 2. Hence the prior art
can be cast as special cases of our more general formulation by imposing
certain limiting restrictions.

###  3.2 Optimization for BiFactor MTL

We propose an efficient learning algorithm for solving the above objective
functionBiFactor MTL. Consider an alternating minimization algorithm, where we
learn the shared representation $\mathbf{F}$ while fixing the task structure
$\mathbf{G}$ and we learn the task structure $\mathbf{G}$ while fixing the
shared representation $\mathbf{F}$. We repeat these steps until we converge to
the locally optimal solution.

Optimizing w.r.t $\mathbf{F}$ gives an equation called generalized Sylvester
equation of the form $AQB^{\top}+CQD^{\top}=E$ for the unknown $Q$. We will
show in the next section on how to solve these linear equation efficiently.
From the objective function, we have:

$\sum_{t}(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})\mathbf{F}(G_{t}^{\top}G_{t})+\lambda_{1}\mathbf{\Sigma}^{-1}F=\sum_{t}\mathbf{X}_{t}^{\top}Y_{t}G_{t}$
(8)

Optimizing w.r.t $\mathbf{G}$ for squared error loss results in the similar
linear equation:

$(\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}\mathbf{X}_{t}\mathbf{F})G_{t}+\lambda_{1}\mathbf{\Omega}^{-1}\mathbf{G}=\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}Y_{t}$
(9)

Optimizing w.r.t $\mathbf{\Omega}$ and $\mathbf{\Sigma}$: The optimization of
the above function w.r.t $\mathbf{\Omega}$ and $\mathbf{\Sigma}$ while fixing
the other unknowns can be learned easily with the following closed-form
solutions (Zhang & Yeung, 2014):

$\displaystyle\mathbf{\Omega}=\frac{(\mathbf{G}\mathbf{G}^{\top})^{\frac{1}{2}}}{tr((\mathbf{G}\mathbf{G}^{\top})^{\frac{1}{2}})}$

$\displaystyle\mathbf{\Sigma}=\frac{(\mathbf{F}\mathbf{F}^{\top})^{\frac{1}{2}}}{tr((\mathbf{F}\mathbf{F}^{\top})^{\frac{1}{2}})}$

###  3.3 TriFactor MTL

As mentioned earlier, one of the restrictions in BiFactor MTL and factored
models is that both the number of feature clusters and task clusters should be
set to $K$. This poses a serious model restriction, by assuming both the
latent task and feature representation live in a same subspace. Such
assumption can significantly hinder the flexibility of the model search space
and we address this problem with a modification to our previous framework.

Following the previous work in matrix tri-factorization, we introduce an
additional factor $\mathbf{S}$ such that we write $\mathbf{W}$ as
$\mathbf{F}\mathbf{S}\mathbf{G}^{\top}$ where $\mathbf{F}$ is a feature
cluster matrix of size $P\times K_{1}$ with $K_{1}$ feature clusters and
$\mathbf{G}$ is a task cluster matrix of size $T\times K_{2}$ with $K_{2}$
task clusters and $\mathbf{S}$ is the matrix that maps feature clusters to
task clusters. With this representation, latent features lie in a $K_{1}$
dimensional subspace and the latent tasks lie in a $K_{2}$ dimensional
subspace.

$\displaystyle\operatorname*{arg\,min}_{\tiny\begin{subarray}{c}\mathbf{F}\in\mathcal{R}^{P\times
K_{1}},\\\ \mathbf{G}\in\mathcal{R}^{T\times K_{2}},\\\
\mathbf{S}\in\mathcal{R}^{K_{1}\times K_{2}}\\\ \mathbf{\Sigma}\succeq
0,\mathbf{\Omega}\succeq 0\end{subarray}}$
$\displaystyle\sum_{t\in[T]}\lVert{Y_{t}-\mathbf{X}_{t}\mathbf{F}\mathbf{S}G_{t}^{\top}}\rVert_{2}^{2}$
(10)
$\displaystyle+\lambda_{1}tr({\mathbf{F}^{\top}\mathbf{\Sigma}^{-1}\mathbf{F}})+\lambda_{2}tr({\mathbf{G}^{\top}\mathbf{\Omega}^{-1}\mathbf{G}})$

The cluster mapping matrix $\mathbf{S}$ introduces an additional degree of
freedom in the factored models and addresses the realistic assumptions
encountered in many applications. Note that we do not consider any
regularization on $\mathbf{S}$ in this paper, but one may impose additional
constraint on $\mathbf{S}$ such as $\ell_{1}$ (sparse penalty), $\ell_{2}^{2}$
(ridge penalty), non-negative constraints, etc, to further improve
performance.

###  3.4 Optimization for TriFactor MTL

We introduce an efficient learning algorithm for solving TriFactor MTL,
similar to the optimization procedure for BiFactor MTL. As before, we consider
an alternating minimization algorithm, where we learn the shared
representation $\mathbf{F}$ while fixing the $\mathbf{G}$ and $\mathbf{S}$, we
learn the task structure $\mathbf{G}$ while fixing the $\mathbf{F}$ and
$\mathbf{S}$ and we learn the cluster mapping matrix $\mathbf{S}$, by fixing
$\mathbf{F}$ and $\mathbf{G}$. We repeat these steps until we converge to a
locally optimal solution.

Optimizing w.r.t $\mathbf{F}$ gives a generalized Sylvester equation as
before.

$\sum_{t}(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})\mathbf{F}(\mathbf{S}G_{t}^{\top}G_{t}\mathbf{S}^{\top})+\lambda_{1}\mathbf{\Sigma}^{-1}\mathbf{F}=\sum_{t}\mathbf{X}_{t}^{\top}Y_{t}G_{t}\mathbf{S}^{\top}$
(11)

Optimizing w.r.t $\mathbf{G}$ gives the following linear equation:

$(\mathbf{S}^{\top}\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}\mathbf{X}_{t}\mathbf{F}\mathbf{S})G_{t}+\lambda_{1}\mathbf{\Omega}^{-1}\mathbf{G}=\mathbf{S}^{\top}\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}Y_{t}$
(12)

for all $t\in[T]$.

Optimizing w.r.t $\mathbf{S}$: Solving for $\mathbf{S}$ results in the
following equation:

$\sum_{t}(\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}\mathbf{X}_{t}\mathbf{F})\mathbf{S}(G_{t}^{\top}G_{t})=\sum_{t}\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}Y_{t}G_{t}$
(13)

Optimizing w.r.t $\mathbf{\Omega}$ and $\mathbf{\Sigma}$: The optimization of
the above function w.r.t $\mathbf{\Omega}$ and $\mathbf{\Sigma}$ while fixing
the other unknowns can be learned as in BiFactorMTL. Note that one may
consider $\ell_{1}$ regularization on $\mathbf{\Omega}$ and $\mathbf{\Sigma}$
to learn the sparse relationship between the tasks and the features (Zhang &
Schneider, 2010).

###  3.5 Solving the Generalized Sylvester Equations

We give some details on how to solve the generalized Sylvester equations
(8,9,11,12,13) encountered in BiFactor and TriFactor MTL optimization steps.
The generalized Sylvester equation of the form $AQB^{\top}+CQD^{\top}=E$ has a
unique solution $Q$ under certain regularity conditions which can be exactly
obtained by an extended version of the classical Bartels-Stewart method whose
complexity is $\mathcal{O}((p+q)^{3})$ for $p\times q$-matrix variable $Q$,
compared to the naive matrix inversion which requires
$\mathcal{O}(p^{3}q^{3})$.

Alternatively one can solve the linear equation using the properties of the
Kronecker product: $(B^{\top}\otimes A)vec(Q)+(D^{\top}\otimes
C)vec(Q)=vec(E)$ where $\otimes$ is the Kronecker product and $vec(.)$
vectorizes $Q$ in a column oriented way. Below, we show the alternative form
for TriFactor MTL equations:

$\displaystyle\sum_{t}((\mathbf{S}G_{t}^{\top}G_{t}\mathbf{S}^{\top})\otimes(\mathbf{X}_{t}^{\top}\mathbf{X}_{t}))vec(\mathbf{F})$
$\displaystyle+\lambda_{1}(\mathbf{I}_{K_{1}}\otimes\mathbf{\Sigma}^{-1})vec(\mathbf{F})$
(14)
$\displaystyle=vec\big{(}\sum_{t}\mathbf{X}_{t}^{\top}Y_{t}G_{t}\mathbf{S}^{\top}\big{)}$

$\displaystyle
diag(\mathbf{S}^{\top}\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}\mathbf{X}_{t}\mathbf{F}\mathbf{S})_{t=1}^{T}vec(\mathbf{G})$
$\displaystyle+(\mathbf{I}_{K_{2}}\otimes\mathbf{\Omega}^{-1})vec(\mathbf{G})$
(15)
$\displaystyle=vec([\mathbf{S}^{\top}\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}Y_{t}]_{t=1}^{T})$

$\displaystyle\sum_{t}((G_{t}^{\top}G_{t})\otimes(\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}\mathbf{X}_{t}\mathbf{F}))vec(S)=vec\big{(}\sum_{t}\mathbf{F}^{\top}\mathbf{X}_{t}^{\top}Y_{t}\big{)}$
(16)

We can do the same for BiFactor MTL, enabling us to use conjugate gradient
descent (CG) to learn our unknown factors whose complexity depends on the
condition number of the matrix $(B^{\top}\otimes A)+(D^{\top}\otimes C)$. To
optimize $\mathbf{F}$, $\mathbf{G}$ and $\mathbf{S}$, we iteratively run
conjugate gradient descent for each factor while fixing the other unknowns
until a convergence condition (tolerance $\leq 10^{-6}$) is met. In addition,
CG can exploit the solution from the previous iteration, low rank structure in
the equation and the fact that the matrix vector products can be computed
relatively efficiently. From our experiments. We find that our algorithm
converges fast, i.e. in a few iterations.

##  4 Experiments

In this section, we report on experiments on both synthetic datasets and three
real world datasets to evaluate the effectiveness of our proposed MTL methods.
We compare both our models with several state-of-the-art baselines discussed
in Section 2. We include the results for Shared Multitask learning (SHAMO)
(Crammer & Mansour, 2012), which uses a K-means like procedure that
simultaneously clusters different tasks using a small pool of $m\ll T$ shared
model. Following (Barzilai & Crammer, 2015), we use gradient-projection
algorithm to optimize the dual of the objective function (Equation 6). In
addition, we compare our results with Single-task learning (STL), which learns
a single model by pooling together the data from all the tasks and Independent
task learning (ITL) which learns each task independently.

Table 1: Performance results (RMSE) on synthetic datasets. The table reports the mean and standard errors over $5$ random runs. The best model and the statistically competitive models (by paired t-test with $\alpha=0.05$) are shown in boldface. Model | syn1 | syn2 | syn3 | syn4 | syn5
---|---|---|---|---|---
STL | 4.79 (0.04) | 5.71 (0.05) | 5.5 (0.04) | 4.02 (0.02) | 5.72 (0.06)
ITL | 1.98 (0.08) | 2.10 (0.09) | 2.01 (0.06) | 1.95 (0.06) | 2.14 (0.07)
SHAMO | 3.63 (0.22) | 4.37 (0.27) | 3.56 (0.23) | 2.76 (0.13) | 4.27 (0.31)
MTFL | 1.91 (0.08) | 1.95 (0.07) | 1.64 (0.06) | 1.47 (0.05) | 1.91 (0.07)
GO-MTL | 1.84 (0.10) | 1.90 (0.08) | 1.72 (0.04) | 1.63 (0.06) | 1.84 (0.06)
BiFactorMTL | 1.85 (0.08) | 1.85 (0.06) | 1.68 (0.08) | 1.37 (0.08) | 1.74 (0.07)
TriFactorMTL | 1.78 (0.08) | 1.83 (0.07) | 1.46 (0.05) | 1.31 (0.02) | 1.68 (0.10)

The parameters for the proposed formulations and several state-of-the-art
baselines are chosen from $3$-fold cross validation. We fix the value of
$\lambda_{1}$ to $0.1$ in order to reduce the search space. The value for
$\lambda_{2}$ is chosen from the search grid
$\\{10^{-3},10^{-2},\ldots,10^{2},10^{3}\\}$. The value for $K$, $K_{1}$ and
$K_{2}$ are chosen from $\\{2,3,5,7,9,10,15\\}$. We evaluate the models using
Root Mean Squared Error (RMSE) for the regression tasks and using $F$-measure
for the classification tasks. For our experiments, we consider the squared
error loss for each task. We repeat all our experiments $5$ times to
compensate for statistical variability. The best model and the statistically
competitive models (by paired t-test with $\alpha=0.05$) are shown in
boldface.

###  4.1 Synthetic Data

We evaluate our models on five synthetic datasets based on the assumptions
considered in both the baselines and the proposed methods. We generate $100$
examples from $\mathbf{X}_{t}\sim\mathcal{N}(0,\mathbf{I}_{P})$ with $P=20$
for each task $t$. All the datasets consist of $30$ tasks with $25$ training
examples per task. Each task is constructed using
$Y_{t}=\mathbf{X}_{t}W_{t}+\mathcal{N}(0,1)$. The task parameters for each
synthetic dataset is generated as follows:

  1. 1.

syn1 dataset consists of $3$ groups of tasks with $10$ tasks in each group
without any overlap. We generate $K=15$ latent features from
$\mathbf{F}\sim\mathcal{N}(0,1)$ and each $W_{t}$ is constructed from linearly
combining $5$ latent features from $\mathbf{F}$.

  2. 2.

syn2 dataset is generated with $3$ overlapping groups of tasks. As before, we
generate $K=15$ latent features from $\mathbf{F}\sim\mathcal{N}(0,1)$ but
tasks in group $1$ are constructed from features $1-7$, tasks in group $2$ are
constructed from features $4-12$ and the tasks in group $3$ are constructed
from features $9-15$.

  3. 3.

syn3 dataset simulates the BiFactor MTL. We randomly generate task covariance
matrix $\mathbf{\Omega}$ and feature covariance matrix $\mathbf{\Sigma}$. We
sample $\mathbf{F}\sim\mathcal{N}(0,\mathbf{\Sigma})$ and
$\mathbf{G}\sim\mathcal{N}(0,\mathbf{\Omega})$ and compute
$\mathbf{W}=\mathbf{F}\mathbf{G}^{\top}$.

  4. 4.

syn4 dataset simulates the TriFactor MTL. We randomly generate task covariance
matrix $\mathbf{\Omega}$ and feature covariance matrix $\mathbf{\Sigma}$. We
sample $\mathbf{F}\sim\mathcal{N}(0,\mathbf{\Sigma})$,
$\mathbf{G}\sim\mathcal{N}(0,\mathbf{\Omega})$ and
$\mathbf{S}\sim\mathcal{U}(0,1)$. We compute the task weight matrix by
$\mathbf{W}=\mathbf{F}\mathbf{S}\mathbf{G}^{\top}$.

  5. 5.

syn5 dataset simulates the experiment with task weight matrix drawn from a
matrix normal distribution (Zhang & Schneider, 2010). We randomly generate
task covariance matrix $\mathbf{\Omega}$ and feature covariance matrix
$\mathbf{\Sigma}$. We sample
$vec(\mathbf{W})\sim\mathcal{N}(0,\mathbf{\Sigma}\otimes\mathbf{\Omega})$.

Table 2: Performance results (RMSE) on school datasets. The table reports the mean and standard errors over $5$ random runs.  Models | STL | ITL | SHAMO | MTFL | GO-MTL | BiFactor | TriFactor
---|---|---|---|---|---|---|---
20% | 12.19 (0.03) | 12.00 (0.04) | 11.91 (0.05) | 11.25 (0.05) | 11.15 (0.05) | 10.68 (0.08) | 10.54 (0.09)
30% | 12.09 (0.07) | 12.01 (0.05) | 10.92 (0.05) | 10.85 (0.02) | 10.53 (0.10) | 10.38 (0.11) | 10.22 (0.08)
40% | 12.00 (0.10) | 11.88 (0.06) | 11.82 (0.06) | 10.61 (0.06) | 10.31 (0.14) | 10.20 (0.13) | 10.12 (0.10)

We compare the proposed methods BiFactor MTL and TriFactor MTL against the
baselines. We can see in Table 1 that BiFactor and TriFactor MTL outperforms
all the baselines in all the synthetic datasets. STL performs the worst since
it combines the data from all the tasks. We can see that the SHAMO performs
better than STL but worse than ITL which shows that learning these tasks
separately is beneficial than combining them to learn a fewer models.

As mentioned earlier, since MTFL is similar to FMTL in Equation 2, we can see
how the results of BiFactor MTL improve when it learns both the task
relationship matrix and the feature relationship matrix. Note that the syn1
and syn2 datasets are based on assumptions in GO-MTL, hence, it performs
better than the other baselines. BiFactor MTL and TriFactor MTL models are
equally competent with GO-MTL which shows that our proposed methods can easily
adapt to these assumptions. Synthetic datasets syn3, syn4 and syn5 are
generated with both the task covariance matrix and the feature covariance
matrix. Since both BiFactor MTL and TriFactor MTL learns task and feature
relationship matrix along with the task weight parameters, they performs
significantly better than other baselines.

###  4.2 Exam Score Prediction

We evaluate the proposed methods on examination score prediction data, a
benchmark dataset in multitask regression reported in several previous
articles (Argyriou et al., 2008; Kumar & Daume, 2012; Zhang & Yeung, 2014)
333http://ttic.uchicago.edu/~argyriou/code/mtl_feat/school_splits.tar. The
school dataset consists of examination scores of $15,362$ students from $139$
schools in London. Each school is considered as a task and we need to predict
exam scores for students from these $139$ schools. The feature set includes
the year of the examination, four school-specific and three student-specific
attributes. We replace each categorical attribute with one binary variable for
each possible attribute value, as suggested in (Argyriou et al., 2008). This
results in $26$ attributes with an additional attribute to account for the
bias term.

Clearly, the dataset has the school and student specific feature clusters that
can help in learning the shared feature representation better than the other
factored baselines. In addition, there must be several task clusters in the
data to account for the differences among the schools. The training and test
sets are obtained by dividing examples of each task into many small datasets,
by varying the size of the training data with $20\%$, $30\%$ and $40\%$, in
order to evaluate the proposed methods on many tasks with limited numbers of
examples.

Table 2 shows the experimental results for school data. All the factorized MTL
methods outperform STL and ITL. We can see that both TriFactor MTL and
BiFactor MTL outperform other baselines significantly. It is interesting to
see that TriFactor MTL performs considerably well even when the tasks have
limited numbers of examples. When there is more training data, the result the
advantage of TriFactor MTL over the strongest baseline GO-MTL is reduced.

###  4.3 Sentiment Analysis

We follow the experimental setup in (Crammer & Mansour, 2012; Barzilai &
Crammer, 2015) and evaluate our algorithm on product reviews from
amazon444http://www.cs.jhu.edu/~mdredze/datasets/sentiment. The dataset
contains product reviews from $14$ domains such as books, dvd, etc. We
consider each domain as a binary classification task. The reviews are stemmed
and stopwords are removed from the review text. We represent each review as a
bag of $5,000$ unigrams/bigrams with TF-IDF scores. We choose $1,000$ reviews
from each domain and each review is associated with a rating from
$\\{1,2,4,5\\}$. The reviews with rating $3$ is not included in this
experiment as such sentiments were ambiguous and therefore cannot be reliably
predicted. h

We ran several experiments on this dataset to test the importance of learning
shared feature representation and co-clustering of tasks and features. In
Experiment i@, we construct $14$ classification tasks with reviews labeled
positive $(+)$ when rating $<3$ and labeled negative $(-)$ when rating $<3$.
We use $240$ training examples for each task and the remaining for test set.
Since all the tasks are essentially same, ITL perform better than all the
other models (with an $F$-measure of $0.749$) by combining data from all the
other tasks. The results for our proposed methods BiFactor MTL ($0.722$) and
TriFactor MTL ($0.733$) are comparable to that of ITL. See supplementary
material for the results of Experiment<EMAIL_ADDRESS>
For Experiment ii@, we split each domain into two equal sets, from which we
create two prediction tasks based on the two different thresholds: whether the
rating for the reviews is $5$ or not and whether the rating for the reviews is
$1$ or not. Obviously, combining all the tasks together will not help in this
setting. Experiments iii@ and iv@ are similar to Experiment ii@, except that
each task is further divided into $2$ or $3$ sub-tasks.

Experiment v@ splits each domain into three equal sets to construct three
prediction tasks based on three different thresholds: whether the rating for
the reviews is $5$ or not, whether the rating for the reviews is $<3$ or not
and whether the rating for the reviews is $1$ or not. This setting captures
the reviews with different levels of sentiments. As before, we build the
dataset for Experiments vi@ and vii@ by further dividing the three prediction
tasks from Experiment v@ into $2$ or $3$ sub-tasks.

The results from our experiments are reported in Table 3. The first four rows
in the table show the number of tasks in each experiment, number of thresholds
considered for the ratings, number of splits constructed from each domain and
the total number of training examples in each task. The general trend is that
factorized models performs significantly better than the other baselines.
Since MTFL, BiFactorMTL and TriFactorMTL learn feature relationship matrix
$\mathbf{\Sigma}$ in addition to the task parameter, they achieve better
results than CMTL, which considers only the task clusters.

We notice that as we increase the number of tasks, the gap between the
performances of TriFactorMTL and BiFactorMTL (and GO-MTL) widens, since the
assumption that the the number of feature and task clusters $K$ should be same
is clearly violated. On the other hand, TriFactorMTL learns with a different
number of feature and task clusters $(K_{1},K_{2})$ and, hence achieves a
better performance than all the other methods considered in these experiments.

Table 3: Performance results (F-measure) for various experiments on sentiment detection. The table reports the mean and standard errors over $5$ random runs.  Data |  ii@ |  iii@ |  iv@ |  v@ |  vi@ |  vii@
---|---|---|---|---|---|---
Tasks | 28 | 56 | 84 | 42 | 86 | 126
| Thresholds
---
(Splits)
2 (2) | 2 (4) | 2 (6) | 3 (3) | 3 (6) | 3 (9)
Train Size | 120 | 60 | 40 | 80 | 40 | 26
STL | 0.429 (0.002) | 0.432 (0.001) | 0.429 (0.002) | 0.400 (0.002) | 0.399 (0.003) | 0.397 (0.001)
ITL | 0.433 (0.001) | 0.440 (0.002) | 0.431 (0.001) | 0.499 (0.001) | 0.486 (0.002) | 0.479 (0.001)
SHAMO | 0.423 (0.002) | 0.437 (0.006) | 0.429 (0.002) | 0.498 (0.006) | 0.460 (0.002) | 0.496 (0.013)
CMTL | 0.557 (0.016) | 0.436 (0.007) | 0.429 (0.004) | 0.508 (0.002) | 0.486 (0.002) | 0.476 (0.002)
MTFL | 0.482 (0.004) | 0.473 (0.002) | 0.432 (0.007) | 0.522 (0.002) | 0.487 (0.003) | 0.481 (0.002)
GO-MTL | 0.582 (0.012) | 0.526 (0.013) | 0.516 (0.007) | 0.587 (0.004) | 0.540 (0.005) | 0.539 (0.008)
BiFactor | 0.611 (0.018) | 0.561 (0.013) | 0.598 (0.002) | 0.643 (0.013) | 0.578 (0.020) | 0.574 (0.052)
TriFactor | 0.627 (0.008) | 0.588 (0.006) | 0.603 (0.012) | 0.655 (0.013) | 0.606 (0.020) | 0.632 (0.029)
Table 4: Performance results (F-measure) on 20Newsgroups dataset. The table reports the mean and standard errors over $5$ random runs.  Models | Task 1 | Task 2 | Task 3 | Task 4 | Task 5
---|---|---|---|---|---
GO-MTL | 0.42 (0.09) | 0.57 (0.06) | 0.42 (0.04) | 0.47 (0.06) | 0.40 (0.03)
BiFactorMTL | 0.42 (0.09) | 0.60 (0.05) | 0.41 (0.04) | 0.49 (0.03) | 0.36 (0.01)
TriFactorMTL | 0.49 (0.03) | 0.63 (0.02) | 0.54 (0.02) | 0.54 (0.02) | 0.51 (0.02)

###  4.4 Transfer Learning

Finally, we evaluate our proposed models on 20Newsgroups dataset for transfer
learning 555http://qwone.com/~jason/20Newsgroups/. The dataset contains
postings from 20 Usenet newsgroups. As before, the postings are stemmed and
the stopwords are removed from the text. We represent each posting as a bag of
$500$ unigrams/bigrams with TF-IDF scores. We construct $10$ tasks from the
postings of the newsgroups. We randomly select a pair of newsgroup classes to
build each one-vs-one classification task. We follow the hold-out experiment
suggested by (Raina et al., 2006) for the transfer learning setup. For each of
the $10$ tasks (target task), we learn $\mathbf{F}$ ($\mathbf{F}$ and
$\mathbf{S}$ in case of TriFactorMTL) from the remaining $9$ tasks (source
tasks).

With $\mathbf{F}$ ($\mathbf{F}$ and $\mathbf{S}$) known from the source tasks,
we select $10\%$ of the data from the target task to learn $G_{target}$. This
experiment shows how well the learned latent feature representation from the
source tasks in a $K$-dimensional subspace ($K_{1}$-dimensional subspace for
TriFactorMTL) adapt to the new task. We evaluate our results on the remaining
data from the target task. We select GO-MTL as our baseline to compare our
results. Since CMTL doesn’t explicitly learn $\mathbf{F}$, we did not include
it in this experiment.

Table 4 shows the results for this experiment. We report the first $5$ tasks
here. See supplementary material for the performance results of all the $10$
tasks. We see that both GO-MTL and BiFactorMTL perform almost the same, since
both of them learn the latent feature representation in a $K$-dimensional
space. As is evident from the table, TriFactorMTL outperforms both GO-MTL and
BiFactorMTL, which shows that learning both the factors $\mathbf{F}$ and
$\mathbf{S}$ improves information transfer from the source tasks to the target
task.

##  5 Conclusions

In this paper, we proposed a novel framework for multitask learning that
factors the task parameters into a shared feature representation and a task
structure to learn from multiple related tasks. We formulated two approaches,
motivated from recent work in multitask latent feature learning. The first
(BiFactor MTL), decomposes the task parameters $\mathbf{W}$ into two low-rank
matrices: latent feature representation $\mathbf{F}$ and task structure
$\mathbf{G}$. As this approach is restrictive on the number of clusters in the
latent feature and task space, we proposed a second method (TriFactor MTL),
which introduces an additional degree of freedom to permit different
clusterings in each. We developed a highly scalable and efficient learning
algorithm using conjugate gradient descent and generalized Sylvester
equations. Extensive empirical analysis on both synthetic and real datasets
show that Trifactor multitask learning outperforms the other state-of-the-art
multitask baselines, thereby demonstrating the effectiveness of the proposed
approach.

## References

  * Amit et al. (2007) Amit, Yonatan, Fink, Michael, Srebro, Nathan, and Ullman, Shimon.  Uncovering shared structures in multiclass classification.  In _Proceedings of the 24th international conference on Machine learning_ , pp. 17–24. ACM, 2007.
  * Argyriou et al. (2008) Argyriou, Andreas, Evgeniou, Theodoros, and Pontil, Massimiliano.  Convex multi-task feature learning.  _Machine Learning_ , 73(3):243–272, 2008.
  * Barzilai & Crammer (2015) Barzilai, Aviad and Crammer, Koby.  Convex multi-task learning by clustering.  In _Proceedings of the 18th International Conference on Artificial Intelligence and Statistics (AISTATS-15)_ , 2015.
  * Crammer & Mansour (2012) Crammer, Koby and Mansour, Yishay.  Learning multiple tasks using shared hypotheses.  In _Advances in Neural Information Processing Systems_ , pp. 1475–1483, 2012.
  * Ding & He (2004) Ding, Chris and He, Xiaofeng.  K-means clustering via principal component analysis.  In _Proceedings of the twenty-first international conference on Machine learning_ , pp. 29. ACM, 2004.
  * Ding et al. (2006) Ding, Chris, Li, Tao, Peng, Wei, and Park, Haesun.  Orthogonal nonnegative matrix t-factorizations for clustering.  In _Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pp. 126–135. ACM, 2006.
  * Dinuzzo et al. (2011) Dinuzzo, Francesco, Ong, Cheng S, Pillonetto, Gianluigi, and Gehler, Peter V.  Learning output kernels with block coordinate descent.  In _Proceedings of the 28th International Conference on Machine Learning (ICML-11)_ , pp. 49–56, 2011.
  * Gu & Zhou (2009) Gu, Quanquan and Zhou, Jie.  Co-clustering on manifolds.  In _Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pp. 359–368. ACM, 2009.
  * Jalali et al. (2010) Jalali, Ali, Sanghavi, Sujay, Ruan, Chao, and Ravikumar, Pradeep K.  A dirty model for multi-task learning.  In _Advances in Neural Information Processing Systems_ , pp. 964–972, 2010.
  * Kumar & Daume (2012) Kumar, Abhishek and Daume, Hal.  Learning task grouping and overlap in multi-task learning.  In _Proceedings of the 29th International Conference on Machine Learning (ICML-12)_ , pp. 1383–1390, 2012.
  * Li & Ding (2006) Li, Tao and Ding, Chris.  The relationships among various nonnegative matrix factorization methods for clustering.  In _Data Mining, 2006. ICDM’06. Sixth International Conference on_ , pp. 362–371. IEEE, 2006.
  * Liu et al. (2009) Liu, Jun, Ji, Shuiwang, and Ye, Jieping.  Multi-task feature learning via efficient l 2, 1-norm minimization.  In _Proceedings of the twenty-fifth conference on uncertainty in artificial intelligence_ , pp. 339–348. AUAI Press, 2009.
  * Maurer et al. (2013) Maurer, Andreas, Pontil, Massimiliano, and Romera-Paredes, Bernardino.  Sparse coding for multitask and transfer learning.  In _ICML (2)_ , pp. 343–351, 2013.
  * Raina et al. (2006) Raina, Rajat, Ng, Andrew Y, and Koller, Daphne.  Constructing informative priors using transfer learning.  In _Proceedings of the 23rd International Conference on Machine Learning_ , pp. 713–720. ACM, 2006.
  * Rothman et al. (2010) Rothman, Adam J, Levina, Elizaveta, and Zhu, Ji.  Sparse multivariate regression with covariance estimation.  _Journal of Computational and Graphical Statistics_ , 19(4):947–962, 2010.
  * Sindhwani et al. (2009) Sindhwani, Vikas, Hu, Jianying, and Mojsilovic, Aleksandra.  Regularized co-clustering with dual supervision.  In _Advances in Neural Information Processing Systems_ , pp. 1505–1512, 2009.
  * Sindhwani et al. (2013) Sindhwani, Vikas, Minh, Ha Quang, and Lozano, Aurélie C.  Scalable matrix-valued kernel learning for high-dimensional nonlinear multivariate regression and granger causality.  In _Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence_ , pp. 586–595. AUAI Press, 2013.
  * Swirszcz & Lozano (2012) Swirszcz, Grzegorz and Lozano, Aurelie C.  Multi-level lasso for sparse multi-task regression.  In _Proceedings of the 29th International Conference on Machine Learning (ICML-12)_ , pp. 361–368, 2012.
  * Wang et al. (2011) Wang, Hua, Nie, Feiping, Huang, Heng, and Makedon, Fillia.  Fast nonnegative matrix tri-factorization for large-scale data co-clustering.  In _Proceedings of the 22nd International Joint Conference on Artificial Intelligence_ , pp. 1553, 2011.
  * Weinberger et al. (2009) Weinberger, Kilian, Dasgupta, Anirban, Langford, John, Smola, Alex, and Attenberg, Josh.  Feature hashing for large scale multitask learning.  In _Proceedings of the 26th Annual International Conference on Machine Learning_ , pp. 1113–1120. ACM, 2009.
  * Xu et al. (2015) Xu, Linli, Huang, Aiqing, Chen, Jianhui, and Chen, Enhong.  Exploiting task-feature co-clusters in multi-task learning.  In _Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence_ , pp. 1931–1937. AAAI Press, 2015.
  * Xue et al. (2007) Xue, Ya, Liao, Xuejun, Carin, Lawrence, and Krishnapuram, Balaji.  Multi-task learning for classification with dirichlet process priors.  _The Journal of Machine Learning Research_ , 8:35–63, 2007\.
  * Zhang et al. (2005) Zhang, Jian, Ghahramani, Zoubin, and Yang, Yiming.  Learning multiple related tasks using latent independent component analysis.  In _Advances in neural information processing systems_ , pp. 1585–1592, 2005.
  * Zhang & Schneider (2010) Zhang, Yi and Schneider, Jeff G.  Learning multiple tasks with a sparse matrix-normal penalty.  In _Advances in Neural Information Processing Systems_ , pp. 2550–2558, 2010.
  * Zhang & Yeung (2014) Zhang, Yu and Yeung, Dit-Yan.  A regularization approach to learning task relationships in multitask learning.  _ACM Transactions on Knowledge Discovery from Data (TKDD)_ , 8(3):12, 2014.
  * Zhong & Kwok (2012) Zhong, Wenliang and Kwok, James T.  Convex multitask learning with flexible task clusters.  In _Proceedings of the 29th International Conference on Machine Learning (ICML-12)_ , pp. 49–56, 2012.

## Sensitivity Analysis

Figure 1 shows the hyper-parameter sensitivity analysis for GO-MTL,
BiFactorMTL and TriFactorMTL. As before, we fix $\lambda_{1}=0.1$. GO-MTL and
BiFactorMTL have two hyper-parameters $\lambda_{2},K$ to tune and TriFactorMTL
have three hyper-parameters $\lambda_{2},K_{1}$ and $K_{2}$ to tune. We can
see from the plots that our proposed models yield stable results even when we
change the $K,K_{1}$ and $K_{2}$. On the other hand, GO-MTL results are
sensitive to the values of $\lambda_{2}$, regularization parameter for sparse
penalty on $\mathbf{G}$.

Figure 1: Top: Sensitivity analysis for the regularization parameter
$\lambda_{2}$ when $\lambda_{1}=0.1$ and $K=2$ (left) and number of clusters
$K_{1}$ and $K_{2}$ when $\lambda_{1}=0.1$ and $\lambda_{2}=0.1$ (right)
calculated for syn5 dataset (RMSE). Middle: Sensitivity analysis for school
dataset (RMSE). Bottom: Sensitivity analysis for sentiment detection
(F-measure).

## Additional Results

Tables 5 and 6 show the complete experimental results for sentiment analysis
and transfer learning experiments.

List of one-vs-one classification tasks used in Table 6

    (Task 1) comp.windows.x
    vs comp.os.ms-windows.misc
    (Task 2) soc.religion.christian
    vs rec.sport.hockey
    (Task 3) misc.forsale
    vs talk.politics.guns
    (Task 4)sci.med
    vs rec.autos
    (Task 5) comp.sys.mac.hardware
    vs talk.politics.misc
    (Task 6) sci.space
    vs alt.atheism
    (Task 7) comp.graphics
    vs comp.sys.ibm.pc.hardware
    (Task 8) talk.politics.mideast
    vs sci.electronics
    (Task 9) rec.motorcycles
    vs talk.religion.misc
    (Task 10) rec.sport.baseball
    vs sci.crypt

Table 5: Performance results (F-measure) for various experiments on sentiment detection. The table reports the mean and standard errors over $5$ random runs. Data |  i@ |  ii@ |  iii@ |  iv@ |  v@ |  vi@ |  vii@
---|---|---|---|---|---|---|---
Tasks | 14 | 28 | 56 | 84 | 42 | 86 | 126
| Thresholds
---
(Splits)
1 (1) | 2 (2) | 2 (4) | 2 (6) | 3 (3) | 3 (6) | 3 (9)
Train Size | 240 | 120 | 60 | 40 | 80 | 40 | 26
STL | 0.749 (0.003) | 0.429 (0.002) | 0.432 (0.001) | 0.429 (0.002) | 0.400 (0.002) | 0.399 (0.003) | 0.397 (0.001)
ITL | 0.713 (0.002) | 0.433 (0.001) | 0.440 (0.002) | 0.431 (0.001) | 0.499 (0.001) | 0.486 (0.002) | 0.479 (0.001)
SHAMO | 0.721 (0.005) | 0.423 (0.002) | 0.437 (0.006) | 0.429 (0.002) | 0.498 (0.006) | 0.460 (0.002) | 0.496 (0.013)
CMTL | 0.713 (0.002) | 0.557 (0.016) | 0.436 (0.007) | 0.429 (0.004) | 0.508 (0.002) | 0.486 (0.002) | 0.476 (0.002)
MTFL | 0.711 (0.002) | 0.482 (0.004) | 0.473 (0.002) | 0.432 (0.007) | 0.522 (0.002) | 0.487 (0.003) | 0.481 (0.002)
GO-MTL | 0.638 (0.006) | 0.582 (0.012) | 0.526 (0.013) | 0.516 (0.007) | 0.587 (0.004) | 0.540 (0.005) | 0.539 (0.008)
BiFactorMTL | 0.722 (0.006) | 0.611 (0.018) | 0.561 (0.013) | 0.598 (0.002) | 0.643 (0.013) | 0.578 (0.020) | 0.574 (0.052)
TriFactorMTL | 0.733 (0.006) | 0.627 (0.008) | 0.588 (0.006) | 0.603 (0.012) | 0.655 (0.013) | 0.606 (0.020) | 0.632 (0.029)
Table 6: Performance results (F-measure) on 20Newsgroups dataset. The table reports the mean and standard errors over $5$ random runs. The best model and the statistically competitive models (by paired t-test with $\alpha=0.05$) are shown in boldface. Models | Task 1 | Task 2 | Task 3 | Task 4 | Task 5 | Task 6 | Task 7 | Task 8 | Task 9 | Task 10
---|---|---|---|---|---|---|---|---|---|---
GO-MTL | 0.42 (0.09) | 0.57 (0.06) | 0.42 (0.04) | 0.47 (0.06) | 0.40 (0.03) | 0.37 (0.02) | 0.35 (0.02) | 0.70 (0.01) | 0.38 (0.00) | 0.42 (0.05)
BiFactorMTL | 0.42 (0.09) | 0.60 (0.05) | 0.41 (0.04) | 0.49 (0.03) | 0.36 (0.01) | 0.42 (0.02) | 0.37 (0.02) | 0.64 (0.02) | 0.38 (0.00) | 0.46 (0.04)
TriFactorMTL | 0.49 (0.03) | 0.63 (0.02) | 0.54 (0.02) | 0.54 (0.02) | 0.51 (0.02) | 0.67 (0.01) | 0.47 (0.02) | 0.66 (0.01) | 0.59 (0.03) | 0.62 (0.01)