Tensors often offer more natural representations of data, e.g., consider video, which consists of obviously correlated images over time. You can turn this into a matrix, but it's just not natural or intuitive (what does a factorization of some matrix-representation of video mean?).

Tensors are trending for several reasons:

• our understanding of multilinear algebra is improving rapidly, specifically in various types of factorizations, which in turn helps us to identify new potential applications (e.g., multiway component analysis)
• software tools are emerging (e.g., Tensorlab) and are being welcomed
• Big Data applications can often be solved using tensors, for example recommender systems, and Big Data itself is hot
• increases in computational power, as some tensor operations can be hefty (this is also one of the major reasons why deep learning is so popular now)

I think your question should be matched with an answer that is equally free flowing and open minded as the question itself. So, here they are my two analogies.

First, unless you're a pure mathematician, you were probably taught univariate probabilities and statistics first. For instance, most likely your first OLS example was probably on a model like this: $$y_i=a+bx_i+e_i$$ Most likely, you went through deriving the estimates through actually minimizing the sum of least squares: $$TSS=\sum_i(y_i-\bar a-\bar b x_i)^2$$ Then you write the FOCs for parameters and get the solution: $$\frac{\partial TTS}{\partial \bar a}=0$$

Then later you're told that there's an easier way of doing this with vector (matrix) notation: $$y=Xb+e$$

and the TTS becomes: $$TTS=(y-X\bar b)'(y-X\bar b)$$

The FOCs are: $$2X'(y-X\bar b)=0$$

And the solution is $$\bar b=(X'X)^{-1}X'y$$

If you're good at linear algebra, you'll stick to the second approach once you've learned it, because it's actually easier than writing down all the sums in the first approach, especially once you get into multivariate statistics.

Hence my analogy is that moving to tensors from matrices is similar to moving from vectors to matrices: if you know tensors some things will look easier this way.

Second, where do the tensors come from? I'm not sure about the whole history of this thing, but I learned them in theoretical mechanics. Certainly, we had a course on tensors, but I didn't understand what was the deal with all these fancy ways to swap indices in that math course. It all started to make sense in the context of studying tension forces.

So, in physics they also start with a simple example of pressure defined as force per unit area, hence: $$F=p\cdot dS$$ This means you can calculate the force vector $F$ by multiplying the pressure $p$ (scalar) by the unit of area $dS$ (normal vector). That is when we have only one infinite plane surface. In this case there's just one perpendicular force. A large balloon would be good example.

However, if you're studying tension inside materials, you are dealing with all possible directions and surfaces. In this case you have forces on any given surface pulling or pushing in all directions, not only perpendicular ones. Some surfaces are torn apart by tangential forces "sideways" etc. So, your equation becomes: $$F=P\cdot dS$$ The force is still a vector $F$ and the surface area is still represented by its normal vector $dS$, but $P$ is a tensor now, not a scalar.

Ok, a scalar and a vector are also tensors :)

Another place where tensors show up naturally is covariance or correlation matrices. Just think of this: how to transform once correlation matrix $C_0$ to another one $C_1$? You realize we can't just do it this way: $$C_\theta(i,j)=C_0(i,j)+ \theta(C_1(i,j)-C_0(i,j)),$$ where $\theta\in[0,1]$ because we need to keep all $C_\theta$ positive semi-definite.

So, we'd have to find the path $\delta C_\theta$ such that $C_1=C_0+\int_\theta\delta C_\theta$, where $\delta C_\theta$ is a small disturbance to a matrix. There are many different paths, and we could search for the shortest ones. That's how we get into Riemannian geometry, manifolds, and... tensors.

UPDATE: what's tensor, anyway?

@amoeba and others got into a lively discussion of the meaning of tensor and whether it's the same as an array. So, I thought an example is in order.

Say, we go to a bazaar to buy groceries, and there are two merchant dudes, $d_1$ and $d_2$. We noticed that if we pay $x_1$ dollars to $d_1$ and $x_2$ dollars to $d_2$ then $d_1$ sells us $y_1=2x_1-x_2$ pounds of apples, and $d_2$ sells us $y_2=-0.5x_1+2x_2$ oranges. For instance, if we pay both 1 dollar, i.e. $x_1=x_2=1$, then we must get 1 pound of apples and 1.5 of oranges.

We can express this relation in the form of a matrix $P$:

 2   -1
-0.5  2 

Then the merchants produce this much apples and oranges if we pay them $x$ dollars: $$y=Px$$

This works exactly like a matrix by vector multiplication.

Now, let's say instead of buying the goods from these merchants separately, we declare that there are two spending bundles we utilize. We either pay both 0.71 dollars, or we pay $d_1$ 0.71 dollars and demand 0.71 dollars from $d_2$ back. Like in the initial case, we go to a bazaar and spend $z_1$ on the bundle one and $z_2$ on the bundle 2.

So, let's look at an example where we spend just $z_1=2$ on bundle 1. In this case, the first merchant gets $x_1=1$ dollars, and the second merchant gets the same $x_2=1$. Hence, we must get the same amounts of produce like in the example above, aren't we?

Maybe, maybe not. You noticed that $P$ matrix is not diagonal. This indicates that for some reason how much one merchant charges for his produce depends also on how much we paid the other merchant. They must get an idea of how much pay them, maybe through roumors? In this case, if we start buying in bundles they'll know for sure how much we pay each of them, because we declare our bundles to the bazaar. In this case, how do we know that the $P$ matrix should stay the same?

Maybe with full information of our payments on the market the pricing formulae would change too! This will change our matrix $P$, and there's no way to say how exactly.

This is where we enter tensors. Essentially, with tensors we say that the calculations do not change when we start trading in bundles instead of directly with each merchant. That's the constraint, that will impose transformation rules on $P$, which we'll call a tensor.

Particularly we may notice that we have an orthonormal basis $\bar d_1,\bar d_2$, where $d_i$ means a payment of 1 dollar to a merchant $i$ and nothing to the other. We may also notice that the bundles also form an orthonormal basis $\bar d_1',\bar d_2'$, which is also a simple rotation of the first basis by 45 degrees counterclockwise. It's also a PC decomposition of the first basis. hence, we are saying that switching to the bundles is simple a change of coordinates, and it should not change the calculations. Note, that this an outside constraint that we imposed on the model. It didn't come from pure math properties of matrices.

Now, our shopping can be expressed as a vector $x=x_1 \bar d_1+x_2\bar d_2$. The vectors are tensors too, btw. The tensor is interesting: it can be represented as $$P=\sum_{ij}p_{ij}\bar d_i\bar d_j$$, and the groceries as $y=y_1 \bar d_1+y_2 \bar d_2$. With groceries $y_i$ means pound of produce from the merchant $i$, not the dollars paid.

Now, when we changed the coordinates to bundles the tensor equation stays the same: $$y=Pz$$

That's nice, but the payment vectors are now in the different basis: $$z=z_1 \bar d_1'+z_2\bar d_2'$$, while we may keep the produce vectors in the old basis $y=y_1 \bar d_1+y_2 \bar d_2$. The tensor changes too:$$P=\sum_{ij}p_{ij}'\bar d_i'\bar d_j'$$. It's easy to derive how the tensor must be transformed, it's going to be $PA$, where the rotation matrix is defined as $\bar d'=A\bar d$. In our case it's the coefficient of the bundle.

We can work out the formulae for tensor transformation, and they'll yield the same result as in the examples with $x_1=x_2=1$ and $z_1=0.71,z_2=0$.

This is not an answer to your question, but an extended comment on the issue that has been raised here in comments by different people, namely: are machine learning "tensors" the same thing as tensors in mathematics?

Now, according to the Cichoki 2014, Era of Big Data Processing: A New Approach via Tensor Networks and Tensor Decompositions, and Cichoki et al. 2014, Tensor Decompositions for Signal Processing Applications,

A higher-order tensor can be interpreted as a multiway array, [...]

A tensor can be thought of as a multi-index numerical array, [...]

Tensors (i.e., multi-way arrays) [...]

So in machine learning / data processing a tensor appears to be simply defined as a multidimensional numerical array. An example of such a 3D tensor would be $1000$ video frames of $640\times 480$ size. A usual $n\times p$ data matrix is an example of a 2D tensor according to this definition.

This is not how tensors are defined in mathematics and physics!

A tensor can be defined as a multidimensional array obeying certain transformation laws under the change of coordinates (see Wikipedia or the first sentence in MathWorld article). A better but equivalent definition (see Wikipedia) says that a tensor on vector space $V$ is an element of $V\otimes\ldots\otimes V^*$. Note that this means that, when represented as multidimensional arrays, tensors are of size $p\times p$ or $p\times p\times p$ etc., where $p$ is the dimensionality of $V$.

All tensors well-known in physics are like that: inertia tensor in mechanics is $3\times 3$, electromagnetic tensor in special relativity is $4\times 4$, Riemann curvature tensor in general relativity is $4\times 4\times 4\times 4$. Curvature and electromagnetic tensors are actually tensor fields, which are sections of tensor bundles (see e.g. here but it gets technical), but all of that is defined over a vector space $V$.

Of course one can construct a tensor product $V\otimes W$ of an $p$-dimensional $V$ and $q$-dimensional $W$ but its elements are usually not called "tensors", as stated e.g. here on Wikipedia:

In principle, one could define a "tensor" simply to be an element of any tensor product. However, the mathematics literature usually reserves the term tensor for an element of a tensor product of a single vector space $V$ and its dual, as above.

One example of a real tensor in statistics would be a covariance matrix. It is $p\times p$ and transforms in a particular way when the coordinate system in the $p$-dimensional feature space $V$ is changed. It is a tensor. But a $n\times p$ data matrix $X$ is not.

But can we at least think of $X$ as an element of tensor product $W\otimes V$, where $W$ is $n$-dimensional and $V$ is $p$-dimensional? For concreteness, let rows in $X$ correspond to people (subjects) and columns to some measurements (features). A change of coordinates in $V$ corresponds to linear transformation of features, and this is done in statistics all the time (think of PCA). But a change of coordinates in $W$ does not seem to correspond to anything meaningful (and I urge anybody who has a counter-example to let me know in the comments). So it does not seem that there is anything gained by considering $X$ as an element of $W\otimes V$.

And indeed, the common notation is to write $X\in\mathbb R^{n\times p}$, where $R^{n\times p}$ is a set of all $n\times p$ matrices (which, by the way, are defined as rectangular arrays of numbers, without any assumed transformation properties).

My conclusion is: (a) machine learning tensors are not math/physics tensors, and (b) it is mostly not useful to see them as elements of tensor products either.

Instead, they are multidimensional generalizations of matrices. Unfortunately, there is no established mathematical term for that, so it seems that this new meaning of "tensor" is now here to stay.

Now I actually agree with most of the content of the other answers. But I'm going to play Devil's advocate on one point. Again, it will be free flowing, so apologies...

Google announced a program called Tensor Flow for deep learning. This made me wonder what was 'tensor' about deep learning, as I couldn't make the connection to the definitions I'd seen.

Deep learning models are all about transformation of elements from one space to another. E.g. if we consider two layers of some network you might write co-ordinate $i$ of a transformed variable $y$ as a nonlinear function of the previous layer, using the fancy summation notation:

$y_i = \sigma(\beta_i^j x_j)$

Now the idea is to chain together a bunch of such transformations in order to arrive at a useful representation of the original co-ordinates. So, for example, after the last transformation of an image a simple logistic regression will produce excellent classification accuracy; whereas on the raw image it would definitely not.

Now, the thing that seems to have been lost from sight is the invariance properties sought in a proper tensor. Particularly when the dimensions of transformed variables may be different from layer to layer. [E.g. some of the stuff I've seen on tensors makes no sense for non square Jacobians - I may be lacking some methods]

What has been retained is the notion of transformations of variables, and that certain representations of a vector may be more useful than others for particular tasks. Analogy being whether it makes more sense to tackle a problem in Cartesian or polar co-ordinates.

EDIT in response to @Aksakal:

The vector can't be perfectly preserved because of the changes in the numbers of coordinates. However, in some sense at least the useful information may be preserved under transformation. For example with PCA we may drop a co-ordinate, so we can't invert the transformation but the dimensionality reduction may be useful nonetheless. If all the successive transformations were invertible, you could map back from the penultimate layer to input space. As it is, I've only seen probabilistic models which enable that (RBMs) by sampling.

Here is a lightly edited (for context) excerpt from Non-Negative Tensor Factorization with Applications to Statistics and Computer Vision, A. Shashua and T. Hazan which gets to the heart of why at least some people are fascinated with tensors.

Any n-dimensional problem can be represented in two dimensional form by concatenating dimensions. Thus for example, the problem of finding a non-negative low rank decomposition of a set of images is a 3-NTF (Non-negative Tensor Factorization), with the images forming the slices of a 3D cube, but can also be represented as an NMF (Non-negative Matrix Factorization) problem by vectorizing the images (images forming columns of a matrix).

There are two reasons why a matrix representation of a collection of images would not be appropriate:

1. Spatial redundancy (pixels, not necessarily neighboring, having similar values) is lost in the vectorization thus we would expect a less efficient factorization, and
2. An NMF decomposition is not unique therefore even if there exists a generative model (of local parts) the NMF would not necessarily move in that direction, which has been verified empirically by Chu, M., Diele, F., Plemmons, R., & Ragni, S. "Optimality, computation and interpretation of nonnegative matrix factorizations" SIAM Journal on Matrix Analysis, 2004. For example, invariant parts on the image set would tend to form ghosts in all the factors and contaminate the sparsity effect. An NTF is almost always unique thus we would expect the NTF scheme to move towards the generative model, and specifically not be influenced by invariant parts.

[EDIT] Just discovered the book by Peter McCullagh, Tensor Methods in Statistics.

Tensors display interest properties in unknown mixture identification in a signal (or an image), especially around the notion of the Canonical Polyadic (CP) tensor decomposition, see for instance Tensors: a Brief Introduction, P. Comon, 2014. The field is known under the name "blind source separation (BSS)":

Tensor decompositions are at the core of many Blind Source Separation (BSS) algorithms, either explicitly or implicitly. In particular, the Canonical Polyadic (CP) tensor decomposition plays a central role in identification of underdetermined mixtures. Despite some similarities, CP and Singular Value Decomposition (SVD) are quite different. More generally, tensors and matrices enjoy different properties, as pointed out in this brief introduction.

Some uniqueness results have been derived for third-order tensors recently: On the uniqueness of the canonical polyadic decomposition of third-order tensors (part 1, part 2), I. Domanov et al., 2013.

Tensor decompositions are nodaways often connected to sparse decompositions, for instance by imposing structure on the decomposition factors (orthogonality, Vandermonde, Hankel), and low rank, to accommodate with non-uniqueness.

With an increasing need for incomplete data analysis and determination of complex measurements from sensors arrays, tensors are increasingly used for matrix completion, latent variable analysis and source separation.

Additional note: apparently, the Canonical Polyadic decomposition is also equivalent to Waring decomposition of a homogeneous polynomial as a sum of powers of linear forms, with applications in system identification (block structured, parallel Wiener-Hammerstein or nonlinear state-space models).