R在大多数人心里的口碑就一个字慢。做转录组分析时处理FPKM矩阵、算α多样性时的反复重采样、拟合SARIMA模型时的滚动调参……明明数据量也就几万行Excel都能拖得动R却在那儿吭哧吭哧烧CPU风扇转得比游戏本还凶。我自己从学生时代到工作被这种“R卡顿”折磨过太多次所以决定把R高性能编程这个方向好好写成一个系列。第一篇讲的是基础习惯比如别滥用全局变量、提前分配向量大小这些这一篇作为二我想直接聊实战里最磨人的几件事瓶颈定位、向量化、并行计算、内存墙以及包管理那些让R变慢的隐形因素。这篇东西适合谁看如果你已经会用R做日常分析但总觉得代码跑得没别人快或者数据一变大就开始盯着进度条发呆那这篇应该能帮你省下好几个通宵。1. 性能分析的起点先测量再优化很多人拿到一段慢代码第一反应是去改写法、上并行、换机器。我的经验是先别急先去搞清楚这段代码的时间到底花在哪。R的慢常常不是均匀分布的而是某一行占了90%的时间其他几千行只是陪跑。如果你不去测量纯粹凭感觉优化很可能折腾了半天把一段本来只要0.1秒的代码优化到0.01秒而真正的瓶颈还在那里纹丝不动。1.1 system.time 和 microbenchmark把“感觉”变成数字R自带的system.time()是最朴素的测速工具。你用system.time(...)包住一个表达式它会返回elapsed时间也就是从开始到结束消耗的秒数。这个函数胜在方便但有个问题单次运行的时间波动很大尤其当你的代码涉及垃圾回收、后台进程或者磁盘读取时跑两次可能差出一倍。所以我一般会用它做粗筛判断某段代码到底是“慢”还是“非常慢”。真正用来做对比的时候我会用microbenchmark包。它可以对同一个操作重复跑很多次自动给出最小值、中位数、平均值还能做统计检验。举个例子我想对比一个循环写法和一个向量化写法x - rnorm(1e6) library(microbenchmark) microbenchmark( loop_way { y - numeric(length(x)) for (i in seq_along(x)) { y[i] - x[i] * 2 1 } }, vec_way x * 2 1, times 50 )在大多数机器上vec_way的中位数可能只有loop_way的几十分之一。这种量级差距一旦被可视化出来你自己都会震惊原来R里的for循环能慢到这种程度。这也是为什么我一直强调性能优化首先要量化不要让模糊的“感觉”帮你做决策。1.2 profvis找到那根“最长的针”system.time能告诉你某段代码整体要多久但如果这段代码跨了几十行、几百行你还需要更细粒度的定位。此时我会用profvis包它的用法极其简单library(profvis) profvis({ # 这里放你的慢代码 for (i in seq_len(10)) { df$new - df$old i } })运行之后RStudio浏览器里会打开一个交互式火焰图横向是时间线纵向是调用栈。你能直观看到每一行代码的耗时占比以及每一行分配了多少内存。这个工具最厉害的地方是能暴露“你以为不重要、实际却吃掉大半时间”的代码行。我自己调试过好几次脚本看着火焰图里某个看似无关紧要的paste0()占了30%耗时那种感觉真的又尴尬又兴奋。1.3 一个反直觉案例t-SNE 栈溢出根本不是性能问题分享一个真实案例。我跑单细胞转录组的t-SNE分析时R报错protect(): protection stack overflow。当时第一反应是数据矩阵太大、内存不够差点准备换服务器。但冷静下来用profvis复现了一遍发现栈溢出根本不发生在t-SNE的核心迭代阶段而是在数据预处理的一个函数里——它递归复制稀疏矩阵时触发了深层递归。这个报错的解决方向有两个。一是用R --max-ppsize500000启动R把保护栈调大二是更根本的思路别再喂高维原矩阵。我当时把上万维的表达矩阵先用PCA压到50维再交给t-SNE不仅不报错速度还快了十倍。这个案例说明什么问题很多所谓的“性能问题”本质上是“定位问题”。你不先测量、不看调用栈、不看上下文可能折腾半天都在治标不治本。2. 从循环思维到向量化思维R提速的第一课如果说性能分析是“找病灶”那向量化就是R里面最常见、也最容易见效的“药”。这一节我想把R的向量化思维掰碎了讲因为它和大部分人脑子里的编程直觉是拧着的。2.1 为什么 for 循环在你脑子里挥之不去大多数人学R之前先学过C、Python或者Java。在这些语言里遍历一个数组天然就是写for循环。于是到了R里面对一列数据脑子里浮现的第一方案依然是“对每一行做什么”。但R的设计哲学是向量化整列数据作为一个对象一次传给底层C代码处理。R的for循环之所以慢核心原因是解释器开销——每执行一次循环体R解释器就要做类型检查、环境查找、变量解析。而向量化调用比如x * 2 1是直接把整个向量交给C语言层面的循环那个循环在编译好的代码里跑速度差一两个数量级非常正常。打个比方for循环像你逐条去翻纸质档案每翻一份还要停下来登记一次向量化像把整摞档案直接扫描进系统机器一次性处理完全部条目。效率当然不在一个量级。2.2 向量化实操FPKM 转 TPM 和分组统计拿转录组分析里最常见的FPKM转TPM来举例。TPM的公式是每个基因的FPKM除以所有基因FPKM之和再乘1e6。新手写法往往是fpkm - c(...) # 假设这是一列FPKM值 n - length(fpkm) tpm - numeric(n) for (i in seq_len(n)) { tpm[i] - fpkm[i] / sum(fpkm) * 1e6 }这个写法有个非常隐蔽的性能杀手sum(fpkm)在循环体内被调用了n次。也就是说每次循环都会重新把整个向量从头到尾加一遍复杂度从O(n)变成了O(n^2)。基因数量如果有两万个这个循环的实际计算量是四亿次加法R要是不卡就怪了。正确的写法是一行向量化tpm - fpkm / sum(fpkm) * 1e6这一行代码做了三件事计算总和一次、每个元素除以总和、再乘以1e6。底层的C代码用一个高效的循环完成全部操作而且全程只遍历了一遍数据。当初我看到自己以前写的循环版本时真想穿越回去抽自己两巴掌。分组统计是另一个典型场景。你想按组求均值写循环的话每个组要subset一次原表频繁复制、频繁扫描慢得离谱。用data.table两三行就能把结果算出来后面我会专门讲它。2.3 提速的正确数据结构data.table 为什么快很多人升级R性能的第一步是从base的data.frame换到data.table。data.frame不是不好而是它有一些非常隐性的性能代价。举个例子df$x - df$x * 2这种操作R内部要复制整个列甚至整个data.frame如果你在一个循环里反复这样改每次都在复制一份大对象。而data.table的更新语法dt[, x : x * 2]是在原对象上直接修改——注意这里没有复制整个表它是引用语义。data.table后面我简称dt的另一个杀手锏是分组聚合。dt[, .(avg mean(x)), by g]这种写法在C层面实现的分组逻辑比base R的aggregate快好几个档次。我自己做单细胞、转录组分析时表达矩阵动辄几万基因、几百样本用data.table去做metadata整理和分组统计几乎没有任何等待感。如果你还在用data.frame处理百万行级别的表格我强烈建议你花一个下午把data.table的语法过一遍这个时间投入的回报率属于R圈里数一数二的。3. 并行计算别让多核闲着也别乱用如果你已经做完向量化和数据结构优化代码还是慢那下一个可以考虑的方向是并行计算。但并行不是银弹用错了反而更慢。3.1 R 的单线程宿命与 parallel 的基本盘R的很多操作是单线程的。哪怕是四核八核的游戏本开一个R进程也往往只用一个核在干活。所以当你有一批相互独立的任务——比如跑100个样本的模型、对100个基因做多重检验、多组重采样计算α多样性指数——并行的收益会非常大。base R的parallel包提供了两个派系。第一个是mclapply它在Unix和macOS上可用原理是fork速度很快因为它直接继承当前进程的内存状态。第二个是parLapply它在所有平台可用包括Windows原理是创建PSOCK集群需要手工启动、停止worker进程还要把变量传到每个worker里去。这个区别很多新手不知道看到网上的mclapply示例就直接抄结果在Windows上报错然后一头雾水。我当初就是从Windows上的parLapply开始接触并行的后来换到Linux服务器上才体验到mclapply的丝滑。3.2 foreach、future 和 BiocParallel三个并行风格怎么选parallel包属于base R能用但写起来不够爽。社区里常见的方案有三个foreach doParallel是经典组合它把循环改造成可并行版本library(foreach) library(doParallel) registerDoParallel(cores 4) res - foreach(i 1:100, .combine c) %dopar% { your_function(i) }future包提供了更统一的抽象。你只需要在开头设置plan(multisession)后面的代码不用大改future_lapply就能并行跑library(future) plan(multisession, workers 4) res - future_lapply(1:100, your_function)未来切换后端时只需要改plan那一行从本机并行切到集群代码主体不用动。生物信息学场景我更推荐BiocParallel。它跟Bioconductor生态咬合得很好像DESeq2、edgeR这些包的内部并行机制都是基于它的。你可以直接用bplapplylibrary(BiocParallel) res - bplapply(1:100, your_function, BPPARAM MulticoreParam(4))注意MulticoreParam在Windows上不能用要换成SnowParam。这也是BiocParallel文档里反复强调的坑。3.3 并行最容易翻车的三个坑对我来说并行最深的三个坑每一个都是血泪换来的。第一个是并行开销。如果单次任务本身只需要0.1秒把它拆给4个worker跑启动worker、序列化传数据、收集结果这些开销可能比省下的时间还多。并行不是越多核越好小任务并行反而更慢。第二个是变量可见性。PSOCK集群里的worker看不到你全局环境里的变量。你必须在调用前显式clusterExport把需要的变量打包发给worker。我经常看到新手写parLapply时调用一个自定义函数报错object not found然后愣在那半天不知道哪里错了。解决方案就是在创建集群后记得导出自定义函数和数据。第三个是随机数种子。并行之后每次运行结果不一样对需要可复现的分析来说非常致命。可以用clusterSetRNGStream设定种子或者用future包时传入seed参数保证每次跑出来的随机序列一致。4. 内存墙大数据跑不动的真正元凶有些代码其实算法没问题但就是跑不动一跑就卡死甚至R直接崩掉。这种情况八成是撞上了内存墙。R在内存使用上有很多反直觉的机制搞懂它们你才知道为什么你的数据明明只有几个G内存却用了几十个G。4.1 R 的 copy-on-modify你无意中复制了多少个亿R有一个非常容易忽视的机制copy-on-modify。你修改一个向量时R不会直接改原对象而是先把整个向量复制一份再对新副本做修改。你可以用tracemem亲眼看一下这种复制x - rnorm(1e7) # 大约80MB tracemem(x) x[1] - 0运行后控制台会输出类似“tracemem[0x0000018... - 0x0000018...”的提示意思是发生了对象复制。这意味着什么一个80MB的向量你只是改了第一个元素R却在背后复制了整个80MB内存瞬间翻了倍。这种机制叠加for循环会非常可怕。如果你在一个循环里反复修改同一个data.frame的同一行或同一列每次修改都可能触发整表复制。循环1000次内存被复制1000次不崩才怪。这也是为什么我一直推荐data.table它用引用语义dt[, x : 1]这种更新不走复制路径。4.2 大矩阵、大表格的落地思路arrow、ff 与分块如果数据大到内存装不下那就得换数据存储和处理思路。第一个思路是选对数据结构。比如一个稀疏的基因表达矩阵你用matrix存每个0都占8字节分分钟几十个G但你用稀疏矩阵比如dgCMatrix只存非零元素内存可能只有原来的十分之一。这也是为什么单细胞分析常用的Seurat底层是稀疏矩阵而不是普通matrix。第二个思路是arrow包。它把数据以列式格式放在磁盘上读取时只加载你需要的列和行而不是整个文件。配合dplyr语法你几乎可以像操作内存数据框一样操作磁盘上的大文件library(arrow) ds - open_dataset(data/parquet_dir) result - ds %% filter(表达量 100) %% select(基因, 样本, 表达量) %% collect()这里collect()只在最后把过滤后的结果拉到内存里前面的filter和select都是在磁盘上完成的。第三个思路是分块处理。把一个大任务拆成多个块每块读入、处理、丢弃这样内存峰值始终可控。我之前处理过一张巨大的定量表就是按基因分块算统计量每块处理完立即rm()再gc()整个过程内存稳稳当当。4.3 监控内存gc()、mem_used 与“系统监测”别等到爆了才想起看内存。RStudio右上角的Environment面板会显示每个对象的大小但整体用量的实时监控我习惯用pryr包的mem_used()library(pryr) mem_used()它会告诉你当前R进程已经占用了多少内存。gc()函数则是显式触发垃圾回收同时打印当前内存使用情况。一个实用习惯是在循环里处理大批量文件时每个迭代末尾显式调用gc()避免垃圾对象一次次累积。Rprofmem()可以做内存剖析逐行显示内存分配情况和数据一样量化之后再优化效率高得多。5. 包管理与环境安装慢、加载慢、跑得慢的三连问最后这部分聊聊包管理。表面上看包的安装和加载好像跟“高性能编程”没什么关系但我觉得它恰恰是让R“变慢”的最大隐形杀手之一。环境没配好你连一个分析都跑不起来更谈不上高性能。5.1 装不上的包权限、源、外部依赖先说最经典的install.packages()在Windows上经常报错“Warning in install.packages: lib C:/Program Files/R/R-4.0.2/library is not writable”。原因很简单你把包装到了系统目录而那个目录没有写权限。解决方法是设置一个个人库路径# 查看当前库路径 .libPaths() # 创建个人库 dir.create(C:/Users/你的用户名/Documents/R/library, recursive TRUE, showWarnings FALSE) .libPaths(c(C:/Users/你的用户名/Documents/R/library, .libPaths()))把这段代码写进.Rprofile文件以后每次启动R都自动生效。从此所有包都装到自己的用户目录下不再碰系统目录权限问题一劳永逸。第二个经典问题是“程序包不存在”。比如你输入install.packages(getoptlong)R提示“不存在叫‘getoptlong’这个名字的程辑包”。这通常不是包真的不存在而是它压根不在CRAN仓库里而是在Bioconductor或GitHub。这时候要用对应的工具# Bioconductor的包 if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(getoptlong) # GitHub的包 remotes::install_github(用户名/仓库名)如果一直用install.packages去找非CRAN包当然找不到。第三个是带外部依赖的包。比如老项目里常碰到的rgdal它依赖GDAL、PROJ这套地理库安装时最容易因为外部库版本不对而卡死。新项目我建议直接用sf替代省去一大串编译烦恼。这虽然是环境问题但它浪费的时间足够你学完整个data.table包所以我把它们放到性能优化里一起说。5.2 让安装飞镜像、并行下载、Rtools 编译链安装慢先检查镜像源。国内用户如果默认连国外CRAN下载速度可以用龟速来形容。在.Rprofile里设置镜像源options(repos c(CRAN https://mirrors.tuna.tsinghua.edu.cn/CRAN/))然后重新启动R之后再install.packages就会快非常多。如果还觉得慢install.packages()支持Ncpus参数可以并行编译install.packages(data.table, Ncpus 4)如果用的是BiocManager同样支持BiocManager::install(DESeq2, Ncpus 4)多核编译对data.table、Rcpp这种带原生代码的包尤其明显四核编译和单核编译的时间差距可以达到三倍以上。这里必须提一下编译工具链。data.table、Rcpp、sf这些包从源码安装时需要编译器。Windows上必须装Rtools而且版本要跟你的R版本严格匹配macOS上要装合适的Command Line Tools。没有工具链的话装包时会提示“compilation failed”。这个坑特别坑看起来像网络问题实际上是编译器缺失折腾半天才发现是环境问题。5.3 加载与重载library 的代价和 devtools CtrlR 的代价很多人写脚本习惯性在顶部library十几个包。这会让R启动和脚本运行前的工作变慢。如果你只用某个包里的一个函数完全可以直接写pkg::function()按需加载而不是把整个包都载入内存。尤其像dplyr、data.table这种本身依赖链很长的大包每多加载一个启动时间就多出一截。开发R包时的体验更明显。你用devtools开发包每次改完代码习惯性按CtrlR——这里说的其实是RStudio中重载整个包的快捷键——它会把包重新加载一遍。如果包里的代码量大依赖多每次load_all()都要好几秒一天下来浪费的时间很可观。我的经验是开发调试阶段尽量只对单独的函数做测试不要频繁重载整个包如果必须重载也要注意把耗时的初始化步骤放在.onLoad里并且做得轻量。事实上这个问题的根源和前面的性能优化是一致的先用数据判断再决定要不要付出那几秒的代价。写到最后说点实在的体会R的“慢”并不全怪R很多时候是我们没按它的脾气来。向量化思维、data.table的引用语义、合理的并行策略、对内存机制的理解这四样吃透之后大部分性能问题都能在你自己手里解决。如果还有跑不动的那基本就是算法复杂度或者硬件到头了该换思路换思路该上服务器上服务器。下一篇我大概率会聊Rcpp和C接口那是把R性能推到极致的最后一块拼图。不过在那之前建议你先把手头的代码重新测一遍指不定一个循环改成向量化今晚就能提前下班。