Embedding as a subgroup
There are several ways you could attempt this, and finding the most efficient method for the types of groups you are interested in might require some experimentation. One way would be to start by finding all (conjugacy classes of) subgroups of $G$ of order dividing $|H|$ and then test them in decreasing order of size for being quotients of $H$.
My inclination would be to try first the simple-minded method of just computing all homomorphisms form $H$ to $G$. I tried that in Magma with a couple fo groups and it worked very quickly. There is a Magma function ${\mathtt {Homomorphisms}}$, which computes the homomorphisms from a finitely presented group to a finite group, up to conjugacy of the image. It finds surjective homomorphisms by default, but there is an option to turn that off. Here is a randomish example with $G$ a simple group of order $20160$ and $H$ a group of order $120$.
I am sure you can do a similar computation in GAP, but I am slightly less familiar with the relevant functions, so I expect someone else could help you with that. Or you could write to the GAP forum mailing list.
> G:=PSL(3,4);
> H:=SmallGroup(120,30);
> H:=FPGroup(H);
> time homs:=Homomorphisms(H,G:Surjective:=false);
Time: 0.220
> [ Order(Kernel(h)) : h in homs ];
[ 120, 60, 60, 20, 60, 30, 30, 30, 30, 30, 30, 30, 30, 12, 12, 15, 15,
15, 15, 15, 15 ]
So in this example, computing the homomorphisms took $0.220$ seconds, and the smallest kernels have order $12$.
Here is another example with the same $G$:
> H:=SmallGroup(120,5);
> H:=FPGroup(H);
> time homs:=Homomorphisms(H,G:Surjective:=false);
Time: 0.040
> [ Order(Kernel(h)) : h in homs ];
[ 120, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 ]
The group $H$ there was ${\rm SL}(2,5)$ and it has found homomorphisms to $A_5$.
In fact I see now that $G$ has no subgroups of order $120$, so perhaps this was a bad choice!
Added later: I have experimented some more, and I see now that the approach above is too naive and works very badly for some types of groups $G$. I tried
> G:=DirectProduct(SmallGroup(100,2),SmallGroup(100,5));
> H:=SmallGroup(100,4);
and the homomorphism computation did not finish in ten minutes. However, the alternative approach of first finding the subgroups of $G$ of order dividing $|H|$ worked fine. Here is how I did that.
> S := [s`subgroup : s in Subgroups(G : OrderDividing := Order(H)) ];
> S := Reverse(S);
> H:=FPGroup(H);
> for s in S do
homs := Homomorphisms(H,s);
if #homs ne 0 then "Found homomorphism, kernel size", #Kernel(homs[1]);
break;
end if; end for;
Found homomorphism, kernel size 25
To complement Derek's answer with an example in GAP, the first three blocks of commands from above look almost identical. Note that H is taken in the polycyclic representation (if one need an fp-group, in GAP the command would be Image(IsomorphismFrGroup(H):
gap> G:=PSL(3,4);
Group([ (3,4,5)(7,9,8)(10,14,18)(11,17,20)(12,15,21)(13,16,19), (1,2,6,7,11,3,
10)(4,14,8,15,16,20,13)(5,18,9,19,21,17,12) ])
gap> H:=SmallGroup(120,30);
<pc group of size 120 with 5 generators>
gap> homs:=AllHomomorphismClasses(H,G);; time;
235
gap> List(homs, h -> Order(Kernel(h)));
[ 120, 60, 60, 15, 15, 15, 30, 30, 30, 60, 30, 30, 15, 15, 15, 30, 30, 30, 12,
12, 20 ]
Next, one could use GAP function AllHomomorphism, but it may be more efficient to use AllHomomorphismClasses which will returns representatives of all homomorphisms from $H$ to $G$ up to $H$-conjugacy:
gap> H:=SmallGroup(120,5);
Group([ (1,2,4,8)(3,6,9,5)(7,12,13,17)(10,14,11,15)(16,20,21,24)
(18,22,19,23), (1,3,7)(2,5,10)(4,9,13)(6,11,8)(12,16,20)(14,18,22)(15,19,23)
(17,21,24) ])
gap> homs:=AllHomomorphismClasses(H,G);; time;
97
gap> List(homs, h -> Order(Kernel(h)));
[ 120, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 ]
The Small Groups Library is included in both GAP and MAGMA. It's crucial that not only the syntax is similar, if not identical, but the numbering is the same:
gap> G:=DirectProduct(SmallGroup(100,2),SmallGroup(100,5));
<pc group of size 10000 with 8 generators>
gap> H:=SmallGroup(100,4);
<pc group of size 100 with 4 generators>
For the last part, one could use LatticeByCyclicExtension to give the upper bound on the order of the subgroup, or use SubgroupsSolvableGroup with several kinds of SizeConsiderFunctions. If I will just use SubgroupsSolvableGroup without any constraints, it takes me 40 seconds to run through all representatives of conjugacy classes of subgroups of $G$ and all classes of homomorphisms to find all possible kernel sizes (time is the CPU time of the last command in milliseconds):
gap> S:=SubgroupsSolvableGroup(G);;
gap> S:=Reversed(S);;
gap> kernelsizes:=[];
[ ]
gap> for s in S do
> homs := AllHomomorphismClasses(H,s);
> for h in homs do
> AddSet( kernelsizes, Size( Kernel( h ) ) );
> od;
> od; time;
39546
gap> kernelsizes;
[ 25, 50, 100 ]
Now an alternative approach which uses LatticeByCyclicExtension and looks only for the subgroups of order dividing 100:
gap> l:=LatticeByCyclicExtension( G,
> function(G) return Size(G) in DivisorsInt((Size(H))); end);
<subgroup lattice of <pc group of size 10000 with 8 generators>, 722 classes,
722 subgroups, restricted under further condition l!.func>
gap> classes:=ConjugacyClassesSubgroups(l);;
gap> S:=List(classes,Representative);;
gap> Set(List(S,Size));
[ 1, 2, 4, 5, 10, 20, 25, 50, 100 ]
gap> S:=Reversed(S);;
gap> for s in S do
> homs := AllHomomorphismClasses(H,s);
> phi := First( homs, h -> Size(Kernel(h)) = 25 );
> if phi <> fail then
> Print("Found homomorphism, kernel size 25 :\n", phi, "\n");
> break;
> fi;
> od;time;
Found homomorphism, kernel size 25 :
Pcgs([ f1, f2, f3, f4 ]) -> [ f3*f5, f3*f6, <identity> of ...,
<identity> of ... ]
31
Thus, we found the first of such homomorphisms immediately. Indeed, we have
gap> phi;
[ f1, f2, f3, f4 ] -> [ f3*f5, f3*f6, <identity> of ..., <identity> of ... ]
gap> Size(Kernel(phi));
25