Seems to be a bit tricky. Here a suggestion which , I think, has to be tested with other examples.

Say your function is

f = gamma + 3 x + a y + z + x y + x y z;

and you know the variables are

var = {x, y, z};

Get the constant term, if any, by setting all variables to 0

In[3]:= const = f /. (# -> 0 & /@ var)    
Out[3]= gamma

Make a new function f1 ( f without constant term )

In[4]:= f1 = f - const    
Out[4]= 3 x + a y + x y + z + x y z

Now for each variable set the other ones to zero and expand the (new) function in a series up to order 1

h[x_] := Normal[Series[f1 /. ( #1 -> 0 &) /@ Complement[var, {x}], {x, 0, 1}]]

Do this for every variable

In[8]:= h /@ var
Out[8]= {3 x, a y, z}

and finally get your linearized function f1L

In[9]:= f1L = const + Plus @@ (h /@ var)
Out[9]= gamma + 3 x + a y + z

Certainely you can combine all these steps in a single function

linF[f_, var_] := Module[{},
  con = f /. (# -> 0 & /@ var);
  f1 = f - con;
  con + Plus @@ (Function[x, Normal[Series[f1, {x, 0, 1}] /. (a_. #1 -> 0 &) /@ 
          Complement[var, {x}]]] /@ var)]

Then for example

In[17]:= linF[123 a + 4 x + beta y + 45 x z - 78 x y z^2, {z, y, x}]
Out[17]= 123 a + 4 x + beta y

Hello Dan,

cool. But

In[25]:= linearize[gamma + 3 x - beta z + 123 x y z]
Out[25]= gamma + 3 x

instead of

In[29]:= linF[gamma + 3 x - beta z + 123 x y z, {z, x, y}]
Out[29]= gamma + 3 x - beta z

Thanks so much, Hans. It is a great code.

Best wishes Omar

I will give a suggestion for a procedure which (attempts) to form a Taylor series for a function of several variables in another thread.