Pyth - 17 bytes

Takes all permutations of the prime factorization, then partitions each one and then products all the partitions, then only retains distinct partitions.

{mS-*Md1s./M.p+1P

Test Suite.

Mathematica, 62 bytes

±1={{}}
±n_:=Union@@(Sort/@Append[n/#]/@±#&/@Most@Divisors@n);

Defines a (recursive) unary operator ± which returns the list of partitions.

Brachylog, 16 bytes

>~l:{1<}a.*?,.=o

This is a function (not a full program) which takes a positive number as input and generates all multiplicative partitions of it. (I also avoided using any prime factorization builtins in this solution, mostly because I wasn't sure they'd help; I might try a more builtin-heavy solution too at some point.)

Try it online! (Extra code has been added around the function here to make it into a full program; if you provide the function shown above to TIO directly, it will run the function but not print its output anywhere, which is kind-of useless as a demonstration.)

This program really disappoints me, because much of it is working around bugs in the Brachylog interpreter and deficiencies in its specification, rather than actually solving the problem; but the interpreter is what it is. (Even with the program like this, the interpreter uses much more memory than it logically should, and crashes due to memory exhaustion, but luckily on small problems it manages to produce the desired output first.) In a hypothetical "perfect version of Brachylog" you could just write ~*.o.:{>1}a,, which would be 4 bytes shorter, but I needed to add extra constraints to help the interpreter out a bit. (I don't really like Brachylog much, and would rather stick to Prolog, but it needed similar hints to make the program work and they're much longer to write. So Brachylog it is.)

Explanation:

As usual, a Brachylog program is a set of constraints; by default, the first constraint constrains the input against an unknown (which I'll call A), the second constraint constrains A against a second unknown B, and so on until we reach the output. Some characters, such as {}, can change this general flow, so I use a different set of letters (e.g. X/Y) to represent unknowns in nested predicates.

>       A is smaller than the input
~l      B has length A
  1<    X is 1, Y is larger
:{1<}a  For each element X of B, it corresponds to an element Y of C
.       C, the output, and D are all identical
*       E is the product of D's elements
?       E, the input, and F are all identical
,       There's no constraint between F and G
.       G, the output, and H are all identical
=       H and I are identical, and need to be evaluated early
o       The output can be produced by sorting I

It's still unclear how the program works, so let's try simplifying the constraints a bit. C, D, G, H, and I are all the same (and equal to the output). E and F are also the same (and equal to the input). So our constraints boil down to this:

• A is the length of B and of the output, and is smaller than the input.
• B consists of all 1s, and isn't particularly useful (it's part of the program simply because in the existing Brachylog interpreter, :{1<}a needs its left argument to have a constrained length, or else the interpreter goes into an infinite loop).
• The output consists entirely of numbers greater than 1 (i.e. greater than the corresponding element of B).
• The product of the elements of the output is equal to the input.
• The output is unchanged by sorting it (i.e. is in sorted order).

Incidentally, I didn't explicitly specify that all the elements of the output are integers, something which might seem to be required; however, Brachylog's constraint solver can't handle non-integers, so it'll conveniently produce only the solutions that involve integers.

Clearly, "the length of the output is smaller than the input" is going to be true whenever the output is a multiplicative partition of the input (because 2^x > x for all nonnegative x, i.e. positive 2^x). So we can disregard that constraint; it's only there to give the Brachylog interpreter a working strategy for evaluating the program. The other constraints (that the output is sorted, that its product is the input, and that its elements are all greater than 1) are the definition of a multiplicative partition, and therefore this function is basically just a direct implementation of the question.

Brachylog, 14 bytes

:{$pp~c:*ao}fd

Maltysen answered this question in 17 bytes of Pyth, so I came up with a 16-byte Brachylog solution that worked by translating the specification of the question to Brachylog. While I was doing that, Dennis wrote a 15-byte Jelly solution. So I had to go down to 14 bytes. This is a function that takes the input as an argument, and returns a list of all partitions (rather than a generator, as with my other solution).

Try it online! (Again, if you provide a Brachylog function to TIO, it runs it but doesn't print the output, so I've added an extra w to make this into a full program by printing the list.)

Unlike my other solution, which didn't make much use of "golfing primitives" but rather stated the problem directly, this one ignores pretty much all the power of Brachylog constraints and does its best Jelly impression instead, writing a chain of constraints for which the left argument is already known (and thus the constraints simply act like Jelly monads rather than full-blown constraints). It therefore uses the same algorithm as @Maltysen's Pyth solution, which is probably the easiest way to solve this using typical golfing primitives. (Interestingly, the "just state the problem" solution in my other answer would be shorter if not for bugs/deficiencies in the Brachylog interpreter, despite its lack of use of golfing primitives. Some day I need to write an "improved Brachylog" in order to get a good solution for this kind of problem; as golfing languages go, Brachylog is actually really verbose.)

The program consists of a generator, and a wrapper around it. First, here's an explanation of the generator:

$pp~c:*ao
$p            Prime factor decomposition of the input
  p           Generate all permutations
   ~c         Generate all inverse concatenations (i.e. partitions)
     :*a      Take the product of each list element in each partition
        o     Sort each partition

This almost solves the problem, but we end up generating many of the partitions many times each. So we need a wrapper to deduplicate the solutions:

:{…}f         Convert generator to list
d             Remove duplicate elements

(Brachylog's habit of requiring a : to specify where a function's second argument starts, even though it can only be a single token long, is really unfortunate from a golfing point of view. There's an obvious improvement to the language right there.)

Python2, 198 191 172 180 bytes

from itertools import*
n=input()
for i in range(2,len(bin(n))):
 for P in combinations_with_replacement(range(2,n),i):
  if reduce(lambda a,b:a*b,P)==n:print(P)
print[(n,),()][n<2]

A full program. This could be improved a lot, so suggestions are deeply welcome!

• @FlipTack saved 9 bytes!
• @orlp saved 19 bytes!

Outputs from the range 1 to 31 (inclusive):

(1,)
(2,)
(3,)
(2, 2), (4,)
(5,)
(2, 3), (6,)
(7,)
(2, 4), (2, 2, 2), (8,)
(3, 3), (9,)
(2, 5), (10,)
(11,)
(2, 6), (3, 4), (2, 2, 3), (12,)
(13,)
(2, 7), (14,)
(3, 5), (15,)
(2, 8), (4, 4), (2, 2, 4), (2, 2, 2, 2), (16,)
(17,)
(2, 9), (3, 6), (2, 3, 3), (18,)
(19,)
(2, 10), (4, 5), (2, 2, 5), (20,)
(3, 7), (21,)
(2, 11), (22,)
(23,)
(2, 12), (3, 8), (4, 6), (2, 2, 6), (2, 3, 4), (2, 2, 2, 3), (24,)
(5, 5), (25,)
(2, 13), (26,)
(3, 9), (3, 3, 3), (27,)
(2, 14), (4, 7), (2, 2, 7), (28,)
(29,)
(2, 15), (3, 10), (5, 6), (2, 3, 5), (30,)
(31,)

Jelly, 15 bytes

ÆfŒ!ŒṖ€;/P€€Ṣ€Q

Try it online!

JavaScript (ES6), 74 67 bytes

f=(n,m=2,a=[])=>n>1?m>n?[]:f(n,m+1,a).concat(f(n/m,m,[...a,m])):[a]

for (var i = 1; i < 31; i++) console.log(JSON.stringify(f(i)));

Directly solves the problem recursively: for each integer m from 2 to n, we take each of the partitions of n / m with a minimum element of m (to avoid duplicate partitions) and append m. (For any m that does not divide n, this gives the empty array, as no arrangement of integers multiplies to a decimal.) We define a base case of the empty array for 1 so as to avoid infinite recursion.

Jelly, 14 13 11 bytes

Ḋx³ŒPQP=¥Ðf

Try it online!

I was fairly sure @Dennis's Jelly solution could be improved. Unfortunately, I didn't manage to beat the Brachylog record, but I did manage to tie it. Update: With @Dennis' help, it's improved now; I guess Jelly takes back the crown.

This program is incredibly inefficient, having O(2ⁿ²) performance (which is why the test case above shows it for input 4). It completes quickly on 4, very slowly on 5, and may not be practically possible to run for larger numbers.

Interestingly, the Brachylog was improved by going from a solution that described the problem (which Brachylog is good at) to a solution that used an algorithm based on factorising the input (which Jelly is good at); meanwhile, the Jelly solution was improved by moving away from its strengths, and going back to a solution that just describes the problem.

Explanation:

Ḋx³ŒPQP=¥Ðf
Ḋ              List of integers from 2 to the input (apparently undocumented)
 x³            Make a number of copies of each that's equal to the input
   ŒP          Take all (possibly noncontiguous) subsequences of that list (!)
     Q         Remove duplicates
         Ðf    Filter, keeping elements where:
      P=         their product is equal to {the original input, by default}
        ¥      Parse preceding two links as a unit

Because the output of Ḋx is sorted, each subsequence must also be sorted, and thus we don't have to sort them individually. So the "the same output in different orders is a duplicate" part of the problem, and the "all values in the output are > 1" part of the problem, get solved by the generation. Apart from that, what we're basically doing here is "find all lists for which P=³", which we do (in an incredibly inefficient way) by generating all the lists in question and then filtering out the incorrect ones.

(Clearly, someone needs to go invent a hybrid of Jelly and Brachylog, plus a really good constraint solver, so that we could write something along the lines of {P=³}~ plus some deduplication code, and solve the program in a much shorter length. That might be some distance away, though.)

Python 3, 65 bytes

f=lambda n,k=2,l=[]:k<=n and f(n/k,k,l+[k])+f(n,k+1,l)or[l][n>1:]

Outputs a list of sorted lists. For example f(20) is [[2, 2, 5], [2, 10], [4, 5], [20]].

Python 3 is used for automatic float division.