JavaScript (ES7), 118 bytes

Returns the coefficients in increasing order.

a=>(i=n=g=p=>n>>i?p.map(Math.round):g(q=[i=t=0],n=!a.map(v=>t+=v**.5*(n>>i++&1||-1)),p.map(x=>q[n++]+=(q[n]=x)*t)))`1`

Try it online!

Method

We don't have any fancy polynomial built-ins, so this implements the same formula as Command Master with a couple of basic map()'s.

Given an array \$a\$ of size \$S\$ and starting with \$p_0=1\$, we compute for each \$n\in[0\dots 2^S[\$:

• \$t=\sum_{0\le i<S}(-1)^{b_i+1}\sqrt a_i\$ where \$b_i\$ is the \$i\$-th bit of \$n\$
• for each \$i\in[0\dots n[\$, starting with \$q_0=0\$:
  • \$q_i\gets q_i+p_i t\$
  • \$q_{i+1}\gets p_i\$
• \$p\gets q\$

Commented

a => (                   // a[] = input array
  i =                    // i = input size (set after the first iteration)
  n =                    // n = output size - 1 (initially zero'ish)
  g = p =>               // g is a recursive function taking an array p[]
  n >> i ?               // if n = 2 ** i:
    p.map(Math.round)    //   return p[] with all values rounded (this
                         //   also coerces the trailing '1' to an integer)
  :                      // else:
    g(                   //   do a recursive call:
      q = [              //     q[] = new output array initialized to [0]
        i =              //     reset i to 0
        t = 0            //     t = sum of square roots, initialized to 0
      ],                 //
      n = !a.map(v =>    //     for each value v in a[]:
        t += v ** .5 * ( //       add ± sqrt(v) to t
          n >> i++ & 1   //       using the i-th bit of n to determine
          || -1          //       the sign
        )                //       (increment i afterwards)
      ),                 //     end of map() / reset n to 0
      p.map(x =>         //     for each value x in p[]:
        q[n++] += (      //       add x * t to q[n]
          q[n] = x       //       and set q[n + 1] to x
        ) * t            //       (increment n afterwards)
      )                  //     end of map()
    )                    //   end of recursive call
)`1`                     // initial call to g with p = ['1']

Octave, 56 bytes

f=@(x)poly((2*(dec2bin(0:2^numel(x)-1)-'0')-1)*sqrt(x)')

Input: a list of integers.

Output: a list of coefficients of a polynomial, whose roots are all possible \$\pm\$ combinations of square roots in input.

Try it online!

Scala + scalaz, 163 bytes

a=>{def P[A]=scalaz.Scalaz.powerset[A]
val b=P(a map math.sqrt)map(_.sum*2)
P(b.map(_-b(0)/2)).groupMap(_.size)(_.product).toSeq.sortBy(_._1)map(_._2.sum.round)}

An expression of type Set[Long] => Seq[Double], but the resulting values are all integers. This is more or less the same approach as Command Master's 05AB1E answer.

Wolfram Language (Mathematica), 26 bytes

MinimalPolynomial@Tr@√#&

Try it online!

Outputs a pure function representing a polynomial.

Extracting the polynomial out of a Root obtained from RootReduce would be a little shorter than using MinimalPolynomial, but such a Root could potentially simplify back to a Sqrt.

Python + SymPy, 53 bytes

-6 from xnor

lambda l:minpoly(sum(map(sqrt,l)))
from sympy import*

Function taking a list of the integers and returning a SymPy minimal polynomial.

Attempt This Online!

Python + NumPy, 81 bytes

lambda l:rint(poly(l**.5@(-1)**(r_[:1<<len(l)]>>c_[:len(l)])))
from numpy import*

Attempt This Online!

How?

Generates what I guess with my Galois theory slowly coming back to me I'm now pretty sure are all the conjugates of the sum of square roots (plus some redundant ones if there are redundant square roots in the sum) by flipping or not flipping all possible combinations of signs. Then multiplying all corresponding linear factors (the poly constructor does that), rounding and hoping for the best. Numerical stability looks atrocious, though the test cases which I stole from @Parcly seem to roughly work.

Charcoal, 59 bytes

ＦＸ²Ｌθ⊞υΣＥθ⎇﹪÷ιＸ²λ²±₂κ₂κＥ⊕Ｌυ﹪%.0fΣＥΦＥＸ²Ｌυ⮌⍘λ²⁼ιΣλΠ∨Ｅ⌕Ａλ1§υν¹

Try it online! Link is to verbose version of code. Explanation: Originally based on @Albert.Lang's answer, but uses the formula subsequently noted in @CommandMaster's answer.

ＦＸ²Ｌθ⊞υΣＥθ⎇﹪÷ιＸ²λ²±₂κ₂κ

Calculate all of the conjugates of the sum of square roots.

Ｅ⊕Ｌυ﹪%.0fΣＥΦＥＸ²Ｌυ⮌⍘λ²⁼ιΣλΠ∨Ｅ⌕Ａλ1§υν¹

For each coefficient, takes all of the ways of choosing that many conjugates, takes the product of each set of conjugates, and takes the sum, then rounds to the nearest integer.

56 bytes using the newer version of Charcoal on ATO:

ＦＸ²Ｌθ⊞υΣＥ₂θ⎇﹪÷ιＸ²λ²±κκＥ⊕Ｌυ﹪%.0fΣＥΦＥＸ²Ｌυ⮌⍘λ²⁼ιΣλΠＥ⌕Ａλ1§υν

Attempt This Online! Link is to verbose version of code. Explanation: As above but the newer version of Charcoal can vectorise SquareRoot over a list and can take the Product of an empty list. I could save a further 4 bytes by not rounding the output (the version of Charcoal on TIO has a bug whereby it incorrectly rounds floating point numbers that are up to 5e-16 less than an integer).

Vyxal, 13 bytes

√:NZΠṠ1v"ƒÞƈ⌊

Try it Online! Takes about 15 seconds for [2, 3, 5], hence the T flag to stop it timing out. Outputs the polynomial backwards, because it saves a byte. Uses a similar approach to Command Master's answer.

Ideally, I'd just be able to use the builtin ∆ṙ, which creates a polynomial from a series of roots. Unfortunately, it's bugged for reasons I don't quite understand, so I have to do it manually with convolutions. This takes a while because sympy isn't particularly good at cancelling out zero expressions (which is also the reason for the ⌊ at the end).

√             # Square roots of each input
   Z          # Zipped with
 :N           # roots negated
    Π         # Get all possible combinations
     Ṡ        # Sum them
      1v"     # Append 1 to each
         ƒÞƈ  # Fold on convolution
            ⌊ # Floor to remove trailing not-quite-zeros

05AB1E, 15 bytes

tæO·¤;-æ.¡g}POò

Attempt This Online!

Computes \$(x-\sqrt{a_1}-\sqrt{a_2}-\cdots)(x+\sqrt{a_1}-\sqrt{a_2}-\cdots)(x-\sqrt{a_1}+\sqrt{a_2}-\cdots)(x+\sqrt{a_1}+\sqrt{a_2}-\cdots)\cdots\$ using Viete's formulas and rounds. Has time complexity \$O(2^{2^n})\$.